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ABSTRACT 

We present a dynamical model of supernova feedback which follows the evolution of pres- 
surised bubbles driven by supemovae in a multi-phase interstellar medium (ISM). The bubbles 
are followed until the point of break-out into the halo, starting from an initial adiabatic phase 
to a radiative phase. We show that a key property which sets the fate of bubbles in the ISM 
is the gas surface density, through the work done by the expansion of bubbles and its role in 
setting the gas scaleheight. The multi -phase description of the ISM is essential, and neglecting 
it leads to order of magnitude differences in the predicted outflow rates. We compare our pre- 
dicted mass loading and outflow velocities to observations of local and high-redshift galaxies 
and find good agreement. With the aim of analysing the dependence of the mass loading of the 
outflow, (3 (i.e. the ratio between the outflow and star formation rates), on galaxy properties, 
we embed our model in the galaxy formation simulation, GALFORM, set in the ACDM frame- 
work. We find that a dependence of (3 solely on the circular velocity, as is widely assumed 
in the literature, is actually a poor description of the outflow rate, as large variations with 
redshift and galaxy properties are obtained. Moreover, we find that below a circular velocity 
of « 80kms^^ the mass loading saturates. A more fundamental relation is that between f3 
and the gas scaleheight of the disk, hg, and the gas fraction, /gas, as /3 oc h^'^f^^^, or the 
gas surface density, Sgas, and the gas fraction, as /3 cx Ti^^^^f^^. We find that using the new 
mass loading model leads to a shallower faint-end slope in the predicted optical and near-IR 
galaxy luminosity functions. 
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1 INTRODUCTION 

An outstanding problem in astrophysics is to understand liow 
galaxies form in dark matter halos. The problem is highly non- 
linear: the stellar mass function of galaxies differs substantially 
from the dark matter halo mass function, with the stellar mass func- 
tion being shallower at the low-mass end and steeper at the high- 
mass end than the halo mass function (see iBaughi2006) . The main 
physical drive r of these differences is thought to be gas cooling 
and feedback ("Larson" 1974'; 'Rees & Ostriker 1977; 'White & Rees' 
1978; Dekel & Silk 1986; White & Frenk 1991; Cole et al. 2000; 
Bower et al .1120061 ; ICroton et al.H2006l) . Feedback from supernovae 
(SNe) and active galactic nuclei (AGN) is thought to suppress star 
formation in low and high stellar mass galaxies, res pectively, low- 
ering the cold baryon fraction in these galaxies (e.g. lpukugita et afl 



11998 ; Mandelbaum et al. 2006; Liu et al. 201Q. 

Observations su ggest that SN-d riven outflows are com- 



ferred outflow rate exceeds the star fo rmation rate (SFR; i MartinI 
ll999l ; lMartir]|2005l ; lBouche et al.ll2012h . suggesting that SN feed- 
back potentially has a large impact on galaxy evolution. The out- 
flow rates inferred from absorption line studies correlate with 
galaxy properties such as SFRs and near-ultraviolet to optical 
colours, indicating that the influence of SN feedback might be dif- 
feren tial with SFR and stellar mass (e.g. lMartirill2005l ; lKomei et al.l 
l2012h . Photometric and kinematic observations of atomic hydro- 
gen shells and holes in the interstellar medium (ISM) of local 
galaxies, in addition to SN remnants observed in X-rays and 
radio, imply that SNe lead to the formation of bubbles within 
the ISM and that the mass carried away is large and able to 
substantiall y change the gas reservo irs of galaxies (e .g. iHeilej 
[l979; Maci eiewski et al. 1996; Pidop rvhora et al.ll2007t) . SN feed- 
back is also thought to be respons ible for the metal enrich- 



mon in galaxies (e .g. Heckman et al. 2000; Shaplev et ; 
Rupkeet al.] |2005| ; Schwartz et al. 2006; Weiner et ; 



talj|2003|; 
aLri2009l; 



ment of the intergalactic medium (e .g. Gnedinll 19981; lAeuirre i 
200 ll; lOppenheimer & Davil2006l; iDubois & Tevssieiil2008l : 



Aguirre et 



Putman et alj|2012l for a recent review) 



Sato etal.ll2009l;lche n et al. 201(il; lRubin et alJl20ld ; lBanerii et al.1 
201 ll ; see lVeilleuxet al.ii2005i for a review). In many cases the in- 



Although the importance of SN feedback is clear from obser- 
vations, the wide range of phenomenological models of SN feed- 
back found in the literature reflect the uncertainty in how this 
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process affects the ISM of galaxies and the intergalactic medium 
(IGM). The key questions are how does the mass loading of winds 
driven by SNe, (5 = A/out /SFR (the ratio between the outflow 
rate, Mout, and the SFR), depend on galaxy properties and what is 
the effect of winds on the evolution of galaxies? 

A common assumption made in galaxy formation modelling 
is that the mass loading (sometimes called the 'mass entrainment' 
of the wind) depends exclusively on the energy input by SNe and 
the circular velocity of the galaxy, which is taken as a proxy for 
the depth of th e gravitational potential well (e.g. IWhite & Reed 
ll978l : IWhite & Frenk 1991) . The specific fom of the dependence 
contains adjustable parameters which are set by requiring that the 
model fits observations, such as the stellar mass function or lu- 
minosity function, etc (e.g . Cole et al. 2000; Springel et al. 2001; 
iBenson et alj|2003l ; ICroton et al.l 120 06). Simple, physically moti- 
vated forms for the explicit dependence of /? on Vdic are based 
on arguments which invoke momentum-driven or energy-driven 
winds, corresponding to dependences of /3 oc v^^}^ and /3 oc 
v7^t, respectively (e.g. [silk 1997; Silk 20 03; Hatl.o n et al. 2003; 
iMurrav et al]|2005l ; IStriiiger et al.ll2012l; see lBensonli20r0l for a re- 
view). 

Hydrodynamic simulations commonly assume constant wind 
velocities, adopting a kinetic feedback scheme in which SNe in- 
ject momentum to neighbouring particles, which are assumed 
to become dyna mically decoupled from the other particles for 
a per iod of time dSpringel & Hemquist' '2003'; 'Scannapieco et al.l 
20061; lOalla Vecchia & Schave 2008; Naravanan et al. 2008; 



Schave et aU20ig) . Alternatively, simple scaling relations between 



the outflow velocity an d the halo circular velocity may be assumed 
(e.g. iDave et al.ll201 ll) . These calcul ations can qualitatively r epro- 
duce the properties of disk galaxies (IScannapi eco et al1l2012l) . The 
wind speed is a free parameter in these sim ulations with values of 
f» 300 - 1000 km/s typically used (see lSchave etai]|20ld for 
an analysis of the impact of changing on the predicted evolution 
of the global density of SFR in a hydrodynamical simulation, and 
IScarmapieco et al]|2012l for a comparison between different simu- 
lations). 

However, such a scheme where the wind speed, «w, is con- 
stant fails to reproduce the stellar mass function, suggesting that 
this parametrisation is too effective in intermediate stellar mass 
galaxies, bu t not efficient enough i n low stellar mas s galaxies 
terain etal.ll200 9; Dave et al. 20 Ij; iBower et al.ll2012h. In addi- 



tion to these problems. iBower et al.l ( 20121) and Guo et al. Il l2013h 



show that simple phenomenological recipes for SN feedback are 
not able to explain t he observed shallow low-mass end of the 
stellar mass function (Drory et al' '2005"; 'Marchesini et alj |2009| ; 
lLi& White! 12009 ; Caputi et al. 2011 ; Bielby et al.. 2012.) . A pos- 
sible explanation for this is that such parametrisations do not accu- 
rately describe the complex process of outflows driven by SNe in 
the interstellar medium and their subsequent propagation through 
the hot halo gas around galaxies. 

ICreasev et al.l ( 1201 3[) analysed the effect of a single SN in the 
ISM by simulating a column through the disk of a galaxy with very 
high mass and spatial resolution. Creasey et al. varied the initial 
conditions in the disk with the aim of covering different gas surface 
densities and gas-to-stellar mass ratios, and found that the mass 
outflow rate depends strongly on the local properties of the ISM, 
su ch as the gas surface density. Similar conclusions were reached 
by Hopkins et al J j2012l) in 4 simulations of individual galaxies in- 
cluding different types of feedback in addition to SN feedback. The 
SN feedback scheme used in Hopkins et al. was not fully resolved 
and hence depends on subgrid modelling of the momentum depo- 



sition of the different types of feedback. Regardless of the details 
of each simulation, both studies point to a breakdown of the clas- 
sical parametrisations used for /3. However, since the simulations 
of both Creasey et al. and Hopkins et al. cover a narrow range of 
environments, the generality of their results is not clear. 

In this paper we implement a fully numerical treatment of 
SN feedback due to bubbles inflated by SNe which expand into 
the ISM. We follow the bubbles during the adiabatic and radiative 
phases assuming spherical symmetry, starting in the star-forming 
regions in the ISM and continuing until the bubble breaks out of 
the galactic disk or is confined. The aims of this paper are (i) to 
study the effect of different physical processes on the expansion of 
bubbles, such as the multi-phase ISM, the gravity from stars and 
dark matter (DM), the temporal changes in the ambient pressure, 
etc., and (ii) to extend previous theoretical work by using the new 
dynamical SN feedback model in the cosmological semi-analytic 
model of galaxy formation, GALFORM. Semi-analytic models have 
the advantage of being able to simulate large cosmological volumes 
containing millions of galaxi es over cosm ic epochs and making 
multiwavelength predictions ( lBauglJl2006l) . This approach makes 
it possible to study a wide enough range of properties and epochs 
to reach robust conclusions about the dependence of /? on galaxy 
properties and to characterise the combination of properties that 
best quantifies the mass outflow rates in galaxies. 

Previous dynamical models of SN feedback in the context of 
cosmological galaxy formation have focused on the evol ution of 
bubbles either i n the ISM or th e h ot halo. For ins tance, iLarsonl 
( 1974) (see alsol Monacoll2004al and lShu et al]|2005l) implemented 
analytic solutions for the evolution of bubbles until their break-out 
from the ISM by neglecting gravity, external pressure and tempo- 
ral changes i n the ambient gas. Bertone et al. (2005, 2007) and 
ISamuietal.1 ([2008') followed the evolution of bubbles in the hot 
halo assuming an ad-hoc mass out flow rate and wind velocity from 
the disk into the halo. iDekel & Silkl d 1986 ) implemented a sim- 
pler model which aimed to estimate the mass ejection rate from 
both the ISM and the halo, using analytic solutions for the evo- 
lution of bubbles in the ISM to calculate an average rate of mass 
injection from the ISM into the halo. Efstathiou ( 2000) went a step 
further, implementing bubble evolution in a multi-phase ISM with 
the hot phase dominating the filling factor, using analytic solutions 
for the evolution of adiabatic bubbles. We improve upon previous 
calculations by including the effects of gravity, radiative losses, ex- 
ternal pressure from the diffuse medium and temporal changes in 
the ambient gas on the expansion of bubbles, all embedded in a 
multi-phase medium. We use the information about the radial pro- 
files of galaxies to calculate mass outflow rates locally. In addition 
to the sophistication of our calculation, another key difference in 
our work is that bubbles expand into the warm component of the 
ISM instead of the hot component, as is assumed in some pre- 
vious work. This is motivated by the results from detailed sim- 
ulations and observations in our Galaxy which point to a rather 
small volume filling factor of hot gas, < 20%, with little mass con- 



tained in this gas phase (e.g.lMac Low et 'Zlll989l;lFeiri5ill200ll; 



Ide AviUez & Breitschwerdal2004 see iHaffner e7Zll2009l for a re- 
view on the warm phase of the ISM). 

In this paper we focus on the ejection of gas from the disk and 
do not attempt to model the expansion of bubbles in the hot halo or 
the rate of gas ejection from the halo into the IGM. In paper II, we 
will implement a full model of the expansion of bubbles in the hot 
halo, following a similar approach to that adopted in this paper, and 
analyse the rate at which mass and metals escape the halo and go 
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into the IGM, and how this depends on galaxy and halo properties 
(Lagos, Baugh & Lacey, 2013, in prep.). 

This paper is organised as follows. ij2]describes the dynamical 
model of SN feedback and the evolution of individual bubbles in 
the ISM. i]2.2l describes the calculation the properties of the diffuse 
medium and how we locate giant molecular clouds (GMCs) in the 
disk. In Sj3]we describe how we include the full dynamical model 
of SN feedback in the galaxy formation simulation GALFORM. In 
21 we analyse the properties of bubbles and the mass and metal 
outflow rate, and their dependence on galaxy properties. We also 
present analytic derivations of some of the relations found in this 
work, giving insight into the physics which sets the outflow rate. 
We study the physical regimes of SN feedback and compare with 
observations of mass outflow rates and velocities in galaxies. In 
^we present a new parametrisation of the outflow rate that accu- 
rately describes the full dynamical calculation of SN feedback and 
compare this to parametrisations that are widely used in the liter- 
ature. In ij6]we show how the new SN feedback model affects the 
galaxy luminosity function and the SFR density evolution. We dis- 
cuss our results and present our conclusions in SjT] In Appendix IaI 
we describe how we calculate the recycled fraction and yield from 
supemovae, in Appendix iBl we explain how we calculate the stellar 
and DM mass enclosed by bubbles, and in Appendix|C]we describe 
how we calculate the overall rates of break-out and confinement of 
bubbles in the ISM. 



2 MODELLING SUPERBUBBLE EXPANSION DRIVEN 
BY SUPERNOVAE 

In this section we describe the physical treatment we apply to bub- 
bles and their expansion in the ISM. We consider that galaxies have 
an ISM which is initially characterised by two gas phases: the dif- 
fuse, atomic phase and the dense, molecular phase. The molecular 
gas is assumed to be locked up in GMCs and stars are allowed 
to form only in these regions. We us e the empirical relation pro- 
posed by iBlitz & Rosolowskvl ^20061) which connects the atomic- 
to-molecular surface density ratio to the hydrostatic gas pressure 
(see §[3T|for details). We use the observed molecular star forma- 
tion (SF) rate coefficient, usf, to calculate the rate at which stars 
form from molecular gas (e.g. Bigiel et al. 2008, 2010). 

The onset of star formation in GMCs results in SNe. SNe 
inject mechanical energy and momentum into the surrounding 
medium, which pressurises the immediate region inflating a cav- 
ity of hot gas, called a SNe driven bubble. We follow the evolution 
of the bubbles from an initial adiabatic phase to a possible radiative 
phase. The interiors of bubbles correspond to a third phase in the 
ISM of galaxies: a hot, low density gas phase. Bubbles start their 
expansion conserving energy, but soon after the expansion starts 
(after a cooling time), the interiors of bubbles become radiative. 
Bubbles then enter into a pressure-driven phase, in which the in- 
terior gas is still hot and highly pressurised. Once this interior gas 
cools radiatively, bubbles continue their evolution conserving mo- 
mentum. 

The main considerations we take into account when following 
the evolution of bubbles are: 

• The injection of energy by SNe lasts for a finite period of time, 
which coiTesponds to the lifetime of a GMC. 

• The gravity of stars and dark matter is included and can decel- 
erate the expansion of bubbles. 

• Temporal changes are followed in the atomic, molecular, stel- 



lar and dark matter contents, with bubbles evolving in this dynam- 
ical environment. 

• We allow bubbles to be offset from the centre of the galaxy but 
they are centered on the midplane of the disk. We therefore consider 
local properties when calculating the expansion of bubbles. 

• Metal enrichment in the ISM due to massive stars takes place 
through bubbles. 

• We follow the radiative cooling in the interior of bubbles to 
make an accurate estimate of the transition between the adiabatic 
and radiative stages of bubble evolution. 

We solve the equations describing the evolution of bubbles numer- 
ically to prevent having to apply restrictive assumptions to features 
we would like to test, such as the effect of ambient pressure and 
gravity on the expansion of bubbles. We make three key assump- 
tions when solving for the evolution of bubbles: 

• Star formation taking place in a single GMC gives rise to a 
new generation of SNe. We assume that the group of SNe in a single 
GMC inflate a single bubble. Thus, each bubble is accelerated by 
a number of SNe, the value of which depends on the SFR in the 
GMC and the initial mass function of stars (IMF). 

• We assume bubbles are spherically symmetric. Observations 
of SNe remnants show that the geometry of bubbles is close to 
spherical in most cases (e.g. lGreenll2009l) . This assumption does 
not restrict the level of accuracy that can be added into the equa- 
tions of momentum and energy describing the evolution of bubbles. 

• We assume that bubbles expand only through the diffuse 
atomic medium and that the gas in GMCs is not affected by 
these expanding bubbles. This is motivated by the fact that GMCs 
are characterised by large gas densities which tend to reflect 
the energy carried out by bubbles rathe r than absor bing it (e.g. 
iMcKee & CowidI 1975l ; lElmegreerj| 199 91). In addition. IWalch et al.l 
1 2OI2I) show that at the moment of explosion of massive stars, the 
surrounding gas has al ready been photo-io nised by the radiation 
emitted by those stars. iHopkins et al.l ( l20I2i) show that this effect 
is also present in their simulations of individual galaxies. This im- 
plies that SNe can efficiently accelerate the suiTounding diffuse gas, 
causing the adiabatic expansion of a bubble to last for longer. 

In §[2jT]we describe the three evolutionary stages for a single 
bubble outlined above and give the equations we use to determine 
the mass, radius, velocity and temperature of the expanding bub- 
bles. In ii l2.2l we describe how we estimate the properties of GMCs 
and the diffuse medium, and how we connect these to the global 
properties of galaxies. 



2.1 Expansion of a single bubble 

Let us consider a bubble located at a distance d from the galactic 
centre and expanding in a diffuse medium characterised by density 
pd, velocity dispersion a^, pressure P^, internal energy density Ud 
and metallicity Zg. 

A single GMC has a SFR of tpcMC and lasts for a time 
TiifcGMC- Within the cloud, the rate of SNe events is t^sn^gmc, 
where »7sn is the number of SNe per solar mass of stars formed. 
The latter depends on the IMF adop t ed. I ndividual SNe release 
EsN = 10^^ erg jAmett et al.lll989l : IWoo slev & Weaver 199^). 
With these definitions in mind we set out the equations we use to 
follow the expansion of bubbles in the following three subsections. 



© 2012 RAS, MNRAS 000.[Tll27l 



4 Claudia del P. Lagos et al. 



(i) ad stage 



(2) 



(ii) pds stage 



(4) 




2.1.1 The adiabatic expansion 

The pressure generated by SNe can significantly exceed tliat of tiie 
ISM, producing a hot cavity. When radiative losses are negligible, 
the hot cavity evolves like a stellar wind bubble which cools adia- 
batically. The interior of the bubble is thermalised and its motion 
drives a shock into the ISM and starts to sweep up the surrounding 
gas ( Ostriker & McKee 1988). The inner structure of the bubble 
corresponds to a thick shell of gas swept-up from the ambient in- 
terstellar medium. The top-panel of Fig. [T| shows a schematic of 
the inner structure of bubbles in this stage, which we refer to with 
the label "ad". The internal gas density profile is illustrated in the 
bottom-right comer. 

The bubble at this stage is characterised by kinetic and thermal 
energies Ek and i5th, respectively, a radius R and an expansion 
speed Us = AR/At, which evolve with time. The total mass of 
the bubble, mb, corresponds to the sum of the mass injected by 
SNe, minj, and the swept-up from the diffuse ISM, msw The rate 
of mass injection depends on i/igmc and the fraction of the total 
mass that is returned to the medium by massive stars, 7?sn, through 
rhinj = -RsN'i/'GMC- Explicit expressions for 77SN and iisN are 
given in Appendix lAl 

The expansion of the inflated bubble is described by the equa- 
tions of energy and mass conservation. 




1 (1) I (2) (a) (4) 



P(r) 



(iii) mds stage 



(3) 




Figure 1. Schematic of the inner stracture of bubbles in three of the expan- 
sion stages considered in our dynamical model of SNe (see §|2j- SNe inject 
energy at a rate Ei^^ , at the centre of the bubble and the ambient medium 
sun'ounds the bubble. A schematic of the gas densities as a function of ra- 
dius depicting the inner stracture of the bubble is shown in the bottom right 
of each panel. Top panel: The adiabatic ('ad') stage. The overpressurised re- 
gion initially expands adiabatically, with the density increasing towards the 
edge of the bubble due to the swept-up gas, producing a thick shell. Mid- 
dle panel: The pressure-driven snowplough ('pds') stage. Once the cooling 
time becomes shorter than the expansion time, the internal mass collapses 
to a shell. The interior mass fueled by the injected mass from SNe remains 
adiabatic. The interior accelerates the outer shell through pressure. Bottom 
panel: The momentum-driven snowplough ('mds') stage. Once the cooling 
time in the interior becomes shorter than the expansion time in the 'pds' 
stage, the interior mass collapses to the shell and forms a bubble with a 
cooled, low density interior The mass and energy injected by SNe modify 
directly the motion of the outer shell through momentum injection. 



-r- = -Binj + 47r R-' 



E 
dE 



Ud — Pd 



GMt{R,d) Gmb 



R 



Pt- 



R 



dTTlb 

"dT 



= minj +47rJ? pdVs. 



(1) 
(2) 



(3) 



Here, E is the total energy of the bubble in the adiabatic stage and 
Einj is the energy injection rate from SNe. 

The total stellar plus DM mass enclosed by a bubble is 
Mt{R, d) and the average density of stars and DM within the bub- 
ble is pt. Both terms act to decelerate the expansion of the bubble 
and come from the gravitational term J^'^ p{r) v{r) g{r) dV in the 
energy conservation equation, where Vb is the volume enclosed by 
the bubble. The term Gpt mb /R represents the increase of grav- 
itational energy internal to the bubble due to the expanding shell 
(see Appendix |B] for a description of the calculation of the stellar 
and DM profiles and the mass enclosed in R). Note that here we 
neglect the self-gravity of the bubble, given that mb <S Mt{R, d). 



The ratio E/{m\^v^ 



ke is calculated using a single 



power-law dependence of the velocity and density on the radius 
inside the bubble (p oc r and v oc r), which gives ke =3/4, for a 
rati o of specific heats of 7 = 5/3 (corresponding to a monoatomic 
gas; lOstriker & McKedl 19881) . The energy injection rate is calcu- 
lated from the SNe rate, rysN'i/'GMC, and the mechanical energy 
produced by an individual SN, Esn, 



Esn rjsN i^GMC- 



(4) 



Note that the pressure of the diffuse medium does not affect 
the energy of the bubbles, given that the diffuse ISM is static with 
respect to the bubbles. This means that there is no coherent mo- 
tion in the ISM, only random motions characterised by a velocity 
dispersion ad- 

For the rate of change in the mass internal to the bubble in 
Eq. |3] the right-hand side of the equation coiTesponds to the rate 
at which mass is incorporated from the diffuse medium into the 
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bubble. We also keep track of the swept-up mass, msw, in order to 
subtract it from the diffuse ISM component when solving the SF 
equations (see §[5), 



drriew . „2 
— — ^ 4tt R pd Vs 
at 



(5) 



Metals produced by nucleosynthesis in stars and ejected by 
SNe are added to the hot cavities. The rate of metal injection by 
SNe into the hot cavity depends on the SFR, ipGMC, the SNe metal 
yield, psN, and the metallicity of the gas from which the stars were 
formed, Zg, and is given by mf^j = (psN + i?SN^g)V'GMC- The 
term psnV'gmc corresponds to the newly synthesized metals and 
^^SN^gV'GMC to the metals present in the gas from which stars 
were made (see Appendix IaI for a description of how the recycled 
fraction and yield are calculated). 

The rate of change in the mass of metals in the interior of 
bubbles and in the swept-up gas component are given by: 



dt 



= "linj + ■ 



dt 



(6) 



(7) 



Similarly to Eq.|5] it is possible to isolate the metals that have been 
incorporated into bubbles from the ISM, mf„ . The internal metal- 
licity of a bubble is therefore Zh = m'^/m\y. This way, the enrich- 
ment of the ISM will depend on the rate of bubble confinement and 
break-out. 

The high temperature of the interior of bubbles results in a 
large sound speed, ^ v^, which makes the time for a sound 
wave to cross the interior much shorter than the expansion time. 
This causes the interior to be isobaric, characterised by a mean 
pressure Pb- We calculate the internal bubble pressure, temperature 
(Tb) and cooling time (tcooi), with the latter two properties defined 
just behind the shock at R (see top panel of Fig.[Tll, using 



mm = 

tcool{R) = 



2 

3" 



Eth 



pnin Pb 

KBPh{R) ' 

3fcBTb(P) 
nbA(Tb(P),Zg)' 



(8) 
(9) 

(10) 



Here, the internal pressure of a bubble is calculated from its internal 
energy, u, fca is Boltzmann's constant, the mean molecular weight 
of a fully ionised gas (i.e. internal to the bubble) is p — 0.62, 
m.H is the mass of a hydrogen atom, A(Tb, ^b) corresponds to the 
cooling function and rib = ph{R)/{p ran ) , is the volume number 
density behind the shock. We adopt the cooling function tables of 
[Sutherland & Dopita ( 1993). 

In order to set the correct initial conditions for the expansion 
in the adiabatic phase, we use the an alytic solutions to the set of 
Eas. lII3| given bv lWeaver et al.l ( ll977l) . These analytic solutions are 
obtained by neglecting the pressure and internal energy of the am- 
bient medium, and the gravity of the stellar plus dark matter com- 
ponent and by assuming that the injected mass is small compared to 
the swept-up mass. We do this for an initial short period of time, t', 
which we quantify in terms of the cooling time, t' ^ 0.1 tcooi- At 
t > t',we follow the solution in the adiabatic stage numerically to 
accurately track the transition to the radiative phase. Our results are 
insensitive to the precise values of t', provided that t' < O.Stcooi- 
The properties of bubbles during this early adiabatic period are: 



Rb{t) 



3 

— a 

5 

4tt 



Einj 

Pd 

Pd 



1/5 



,3/5 



1/5 



-2/5 



msw(t) Zg, 

msw(t) + PsNl/'GMC t, 



irih{t) = 

mb(t) = mf„(t) + (psN + PsN^g)V'GMC i, 



(11) 
(12) 

(13) 

(14) 
(15) 
(16) 



where a = 0.86. Egs. 1 151 161 account for the injected metals and 
mass from the dying stars. 



2.1.2 Pressure-driven snowplough expansion 

As the temperature of the bubble decreases with time, the cool- 
ing time becomes sufficiently short so as to be comparable with 
the expansion time of the bubble. At this stage, radiative losses 
from the expanding thick shell can no longer be neglected and 
the shocked swept-up material quickly becomes thermally unsta- 
ble and collapses into a thin, dense shell. The shocked mass ejected 
by SNe in the interior of the thin shell still conserves its energy and 
the bubble enters a pressure-driven phase. The energy injected by 
SNe modifies the thermal energy of the shocked interior We refer 
to properties of bubbles in this stage with the label "pds", denoting 
pressure-driven snowplough (see middle panel of Fig. [T). 

In this phase bubbles are characterised by the swept-up mass 
accumulated in a thin shell, msh, and an interior mass, mint . The 
interior of the bubble is still isobaric, characterised by a mean pres- 
sure. Pint- We consider that the density of the shocked SNe injected 
material is constant and is calculated Pint — ^int /(4/37rP^). 

We calculate Put using Eq.[8] Pint = i^int /2-kR? , where Bint 
is the interior energy of the bubble and is calculated from the energy 
gained from SNe (iJinj) and the energy loss due to the work done 
by the interior gas on the expanding shell. 



dt 



Einj - 47r P" V, Pin 



(17) 



The rate of change of mass and metals in the interior of bubbles are 
set by the mass and metals injection rates by SNe, mint ~ minj 
andm?t = m?j. 

The temperature and cooling time in the interior of the bub- 
ble are calculated following Eqs.|9]and[T0] but replacing p{R) by 
Pint = mint/{|7rP^), Pb by Pnt and Zh by Zint = mfnt/rriint. 

The equations of motion and of the conservation of the total 
mass and mass in metals for the shell in the pressure-driven stage 
are 



dt 

drrish 
dt 



dt 



47r R^ (Pin 



47r R Pd Vs 



— Atv R Pd Vs Zg 



GMt{R,d) 

P2 



nish (18) 
(19) 
(20) 



Note that the expansion of the bubbles is driven by the pressure dif- 
ference (Pint — Pd)- The gravitational term G Mtrriaii/ R? comes 
from integrating g5M over all the mass elements inside a radius 
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that is comoving with the diffuse medium in the equation of mo- 
tion for an element of fluid of mass SM. We neglect the shell self- 
gravity, given that rus <€. Mt{R, d). 



2.2 Properties of molecular clouds and the diffuse medium in 
galaxies 

In this section, we describe how we calculate the properties of 
GMCs and the diffuse medium, and explain the techniques used 
to follow their evolution throughout the ISM. 



2.1.3 Momentum-driven snowplough expansion 

When the expansion time in the pds stage becomes longer than the 
cooling time of the interior, the bubble enters to the momentum- 
driven phase. The cavity interior to the bubble is composed of low- 
density cooled gas of total mass mint. This interior mass corre- 
sponds to the ejected mass from SNe that has not yet had enough 
time to encounter the shell. The explosions at the centre inject mass 
and momentum into the shell. The interior density is calculated 
from the continuity equation. 



Pint = 



(21) 



The density of the ejected material drops with radius and by the 
time the ejected gas encounters the shell, most of the energy input 
by SNe has become kinetic energy. Therefore, SNe ejected mate- 
rial acts on the shell by increasing the momentum of the shell (see 
schematic in the bottom panel of Fig. [TJ. We therefore consider 



that Vinj = Y 2 Einj/rhinj. The equations describing the change of 
mass and mass in metals of the bubble interior are: 



dt 
dt 



(22) 
(23) 



Here, the amount of injected mass that remains in the interior of the 
bubble depends on the velocity ratio Hs/uinj, which means that if 
the shell expands slowly, most of the mass injected by SNe quickly 
reaches the shell. Note that gravity is neglected in the motion of the 
interior material. 

The equations describing the conservation of momentum, total 
mass and mass in metals for the mds stage are. 



dt 



dmsh 
dt 

dm,\ 
dt 



. , , GMt{R,d) 
mini ("inj - "sj ^ msh 



i?2 



-47r iJ^Pd, 



mini 1 1 — ] +471 pd «s 



+4TrR pd^s^d- 



(24) 
(25) 

(26) 



Note that the expansion of the bubbles is driven by the velocity 
gradient (uinj ~ v^). 

If the bubble has a radius which exceeds the scale height of 
the galaxy, part of the bubble would be expanding in a lower den- 
sity medium (see bottom panel of Fig. |2](. We account for this by 
including a correction factor in the density of the diffuse medium 
when _R > /ig, Pd = Pd (1 — hg/R), which accounts for the frac- 
tion of the surface of the bubble outside the disk. We replace pd by 
Pd in the set of equations describing the evolution of bubbles. 



2.2.1 Molecular cloud properties 

The dynamical evolution described above corresponds to a single 
bubble driven by the SF taking place in one GMC. In order to in- 
corporate this evolution into the galaxy formation context, we con- 
sider GMC formation in the ISM of galaxies and subsequent SF 
in GMCs. For this, it is necessary to define the GMC mass, SF 
efficiency and the timescales for the formation and destruction of 
GMCs. We first define individual GMC properties and then con- 
nect them to galaxy properties to estimate their number and radial 
distribution in ^ 12.2.31 

GMC mass. Motivated by observations of the Milky Way 
and nearby galaxies, we consider GMCs t o have typical masses 
of mcMC » 10^ - 10^M(7) (e.g. jSolom on et al.' 'l987|; 
Iwilliams & McKed 1 19971 : I Oka et alj I2OO1I : iRosolowskv & Bliti 
[2005 ). We assume that GMCs are fully molecular and that all the 
molecular gas in galaxies is locked up in GMCs. This is a good 
approximation for most local galaxies, in which more t han 90% 
of the molecular gas is in gravitationally bound clouds j Ferrierd 
2001). However, it is important to note that in the densest nearby 
starburst galaxies, so me molecular gas is also fou nd in the diffuse 
component (e.g. M64: lRosolowskv & BIitzl2005l) . 

The SFR per GMC. V'gmc depends on the GMC mass and 
the molecular SF coefficient rate, i^sF, as "^gmc = J^Spn^GMC- 
To ensure consistency with the global SF law, we use the same SF 
rate coefficient defined in § |2] This implies that, as we incorpo- 
rate the dynamical SNe feedback model in the galaxy formation 
simulation, GMCs forming stars in the disk have different deple- 
tion timescales than those forming stars in the bulge (see 5* 13.11 
for details). This difference in the SF timescales of GMCs in nor- 
mal star-formin g galaxies and starburs ts (SBs) has been proposed 
theoretically bv Krumholz et al] ( l2009l) . Krumholz et al. argue that 
in normal galaxies the ambient pressure is negligible compared to 
the internal pressure of GMCs, and therefore, the properties setting 
the SF are close to universal. However, in high gas density envi- 
ronments appropriate to SBs, the ambient pressure becomes equal 
to the typical GMC pressure, and therefore, in order to maintain 
GMCs as bound objects, their properties need to change accord- 
ing to the ambient pressure. This naturally produces a dichotomy 
between normal star-forming galaxies and starburst galaxies. 

GMC lifetime. The formation and destruction timescales 
of GMCs depend on the properties of the ISM: gas density, 
convergence flow veloci ties, magnetic fields, turbulence, etc. 
jMcKee & Ostriked lIoOTh . GMCs can form through large-scale 
self-gravitating instabilities, which can include Parker, Jeans, 



magneto-Jeans and/o r magneto-rotational instabilities (e.g. Chiezd 
'1987'; Malonev 1988:'Elme greerilll989l : lM"cKee & HollimarJll99g| : 



Krumholz & McKee 2005 ), or through collisions of large-scale 
gas flows (e.g.|B allesteros- Paredes et al.lll999l : iHeitsch et al.ll200^ : 
IVazquez-Semadeni et al.ii2006l) . GMCs in these formation scenar- 
ios tend to last ~ 1 — 3 crossing times before being destroyed 
by stellar feedback (i.e. proto-stellar and stellar winds, and HII 
regions). Observationally, the lifetime of GMCs is inferred from 
statistical relations between the location of G MCs and young star 
clusters and is in the range 10 - 30 Myr (e.g. IBlitz & Shulll980l: 
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lEngargiola etai]|2003l ; Isiitz et alJlioOTh . We therefore restrict the 
range of the lifetimes of GMCs to riife,GMC = 10 — 30 Myrs. 

2.2.2 Properties of the perx'asive interstellar medium 

We assume that the diffuse pervasive medium in the ISM is fully 
atomic. We define the relevant properties of the diffuse medium 
(see Eas lH26| l as a function of radius for the disk and bulge. 

For the gas surface density profiles of the disk and bulge, we 
assume that both are well described by exponential profiles with 
half-mass radii, rso.d and r5o,b, respectively. This is done for sim- 
plicity. However, it has been shown that the neutral gas (atomic plus 
molecu lar) in n earby spiral galaxies follows an exponential radial 
profile teigiel & Blitz 2012). Davis et al. (2012) found that this is 
also the case in a large percentage of early-type galaxies in the local 
Universe. In interacting galaxies and galaxy mergers, Davis et al. 
show that the gas can have very disturbed kinematics, and in these 
cases our approximation is no longer valid. 

To calculate the HI surface density we follo w iLagos et"ail 
( l2011hl) and use the iBlitz & Rosolowskvl j200d) pressure law 
(§ [3j. We assume this pressure-law also holds in higher gas den- 
sity media, typical of SBs. Hydrodynamic simulations including 
the formation of H2 have shown that, for extreme gas densities, 
the relation between hydrostatic pressure and the Ehj/Shi ra- 
tio deviates from the empirical pr essure law resulting in more H2 
jPelupessv & Papadopouloj|2009l ). If the conclusions of Pelupessy 
et al. are correct, our assumption that the Blitz & Rosolowsky law 
holds for SBs would represent an upper limit for the HI mass. The 
effect of this systematic on the final result of SNe feedback is highly 
non-linear given that having more HI mass makes the expansion of 
bubbles more difficult, but in the case of escape, more outflow mass 
is released from the galaxy. 

We assume that gas motions in the diffuse medium are domi- 
nated by a random component an d we choose the v ertical velocity 
dispersion to be (Td = 10 km dLerov et al.l2008l) . The source of 
the motion of the diffuse ISM is not relevant so long as it gives rise 
to gas dominated by random motions. The assumption of random 
moti ons is consistent with turbulence and thermally drive n motions 
(e.g. IWada et"ai]|2002l : [Schav3l2004l : loobbs et al.ll2oTTI) . We esti- 
mate the gaseous disk scale height, the volume density and thermal 
pressure as a function of radius, /ig(ri), Pd(^i) and Pd(''i), respec- 
tively. The set of equations defining these properties is 



pd(n) 



-I 



Satom (^i) 

2ftg(n) ' 

Pd(n)crd. 



+ 



(27) 

(28) 
(29) 



Here is the velocity dispersion of the stars, and Eatom(''i), 
Eg (n) and E* (ri) are the atomic, total gas (molecular plus atomic) 
and stellar surface densities, respectively, at r^. In Appendix IB II we 
describe the calculation of ctj, and the origin of the expression for 
/ig. The choice of (Jd fixes the internal energy of the diffuse medium 
throughout the disk and bulge, so that u = 3/2 Pd- 

Note that we include the contribution of helium in pd (ri). The 
filling factor of molecular clouds in the ISM is very small, typically 
Fgmc ~ 0.01 (iMcKee & Ost rikej|2007h . so we assume that the 
filling factor of the diffuse gas is _Fd = 1 and therefore we do not 
include it in Eq. |27|29l 

The gas scale height includes the gravitational effect of stars 



through E*(ri). The underlying assumption in Eq. |27]is that the 
galaxy is in vertical equilibrium and that the diffuse medium is 
characterised by a uniform pressureQ. 



2.2.3 Connecting GMCs and galaxy properties 

We follow the evolution of bubbles in rings within the disk and the 
bulge, and assume cylindrical symmetry: all bubbles at a given ra- 
dius Ti from the centre are identical, where i — l..Nr. We estimate 
the number of molecular clouds in the ISM at a given timestep that 
give rise to a new generation of bubbles. If at a timestep t — tj the 
radial profile of molecular mass is Enioi('', tj), the total number of 
GMCs in an annulus of radius ri and width 5r is. 



Ng 



2-X 



MCs,i,j = 



ri+(5r/2 , 
ri-Sr/2 ' 



i{r,tj)rdr 



m-GMC 



(30) 



The rate of GMC formation in the annulus i in a given time tj is 
therefore estimated as. 



Ng 



MC,ncw,i,j 



A^GMCs.i.j 
TlifcGMC 



(31) 



Note that by fixing the SF rate coefficient, i^sf, and the properties 
of GMCs, we are implicitly assuming that all GMCs at a given 
timestep are forming stars. 

We performed tests to choose the value of A'^r to ensure con- 
vergence in the results presented in this work. These tests suggests 
A^r ~ 10. The spatial extent of each ring i depends on the total ex- 
tent of the disk we choose to resolve. We model out to Srso in disk 
radius, so the molecular mass enclosed is > 99.999% of the total. 
This defines the extent of the individual annuli, 5r = Srso/A^r- 



2.2.4 Bubble confinement and break-out 

Confinement. If bubbles are slowed down sufficiently, they are as- 
sumed to mix with the surrounding medium. The condition for 
mixing to take place is obtained by comparing the bubble expan- 
sion velocity to the velocity dispersion of the diffuse component of 
the ISM. Confinement takes place if Us ^ Ud. If this happens, we 
assume instantaneous mixing and add the mass and metals of the 
bubble to the diffuse medium of the ISM. 

Break-out from the ISM. If a bubble reaches the edge of the disk 
or the bulge with an expansion velocity exceeding the sound speed 
of the diffuse ISM, it is assumed to break out from the ISM. The 
edge is defined as a fixed fraction of the gas scale height, /r hg 
(see § 12.2.21 for the definition of gas scale height). The opening 
angle of the wind at the moment it escapes from the galaxy is given 
by ^ = 2 arccos(l//i-), assuming that bubbles are centered at the 
midplane of the disk. A fraction /bo of the mass and metals carried 
away by bubbles will escape from the galaxy. This depends on the 
choice of /r — R/h^ is given by 



/bo 



R 



= i-/rv 



(32) 



^ IShettv & Ostrikej l l2012h use a set of vertically resolved hydrodynamic 
simulations to sh ow that vertical equi l ibrium is reached within a vertical 
crossing time and lKovama & Ostrikej j2009l) show that vaiiations in pres- 
sure vertically are within a factor of 2. 
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Table 1. List of parameters in the dynamical SNe feedback model. In the right-hand column, theoretical and observational constraints on these parameters are 
described. The values adopted in our preferred model (referred to as the standard model in the text) are indicated in parentheses. 



symbol 


parameter 


range and value adopted 


constraints from obs. and theory 


GMC pai'ameters 


A^GMC 
nife.GMC 


typical mass of a GMC 
Lifetime of a GMC 


Mgmc = 10^ - 106 Mq 
(std. model A/gmc = 10^ Mq) 

tlifcGMC = 10-30 Myr 
(std. model tiifc,GMC = 10 Myr) 


Estimated bv Solomon et al. (1987). 
Williams & McKee (1997). 
Observations and theoretical arguments 
favour values in the range given here 
(e.g. Blitz & Shu 1980: Dobbs et al. 2011). 


Diffuse medium parameters 




velocity dispersion of 
the gas in disks 


(Td ~ 5 — 10 km s ^ 
(std. model cr^ = lOkms^^) 


Lerov et al. (2008). Used to calculate 
Pd : i^d and hg . 


Disk parameters 


/* 

A 


ratio of the scale radius to 
the scale height of the stellar disk 
Defines the radius at which 
bubbles are assumed to have 
escaped the galaxy. 


/* = rs//lstar ~ 7.3 
/r = 1.1-2 

(std. model /r = 1.5) 


Kreael et al. (2002). Used to calculate 
Pext and hg. 

In principle /r is a free parameter. 
However, we set a range within which 
we vary /r as to get a break-out 
mass fraction consistent with previous 
theoretical estimates 
(e.g. MacLow & McCrav 1988; 
Fuiita et al. 2009). 


SF parameters 


Po 

ap 


SFR coefficient 

Pressure normalisation 

Power-law index in 
pressure law 


i^sF = 0.25 - IGyr-i 
(std. model i/gp = 0.5 Gyr~^) 

log(Po/kB[cm-='K]) = 4.19 - 4.54 
(std. model 
log(Po/kB[cm-3K]) =4.54) 
ap = 0.73 - 0.92 
(std. model ap = 0.92) 


Detennines the SFR per unit 
molecular mass EsFR = ''SF ^mol- 
Measured bv e.g. Lerov et al. (2008). 
Shj/Shi = (Pext/Po)"P. Measured 
bv e.g. Wong & Blitz f2002). BUtz & 
Rosolowskv (2006). Lerov et al. (2008). 
Measured (see authors above). 



A fraction (1 — /bo) of the mass and metals carried away by bub- 
bles is assumed to be confined in the ISM. The physical motivation 
for this choice is that the gas expanding along the major axis of 
the disk does not escape and that, in the case of the gas expanding 
perpendicular to the midplane of the disk, Rayleigh-Taylor instabil- 
ities grow at the edge of the ambient gas due to the drastic change 
of density. These instabilities produce fragmentation in the swept- 
up mass and some of this m ateri al is reincorporated in to the galaxy. 
iMacLow & McCravl ( Il988l) and lMac Low et alj ^1989), by means 
of hydrodynamic al simulations, estimate d /r « 1 — 2 for a Milky 
Way-like galaxy. iMac Low et alj ( 1 19891) show that approximately 
10% of the mass contained in shells at the point of break-out ac- 
celerates upwards and ~ 90% stays in the ISM. Similar values 
have been obtained by more sophisticated hydrodyn amical simu- 
lations (e.g. lde Avillez & Berrvll200lHFuiita et al.ll2009i) . In detail, 
the break-out radius and the mass in shells escaping the galaxy disk 
is thought to mainly depend on the density contrast between the 
disk and halo gas which sets the development of instabilities which 
fragments the bubble shells. Other hydrodynamical effects, such 
as weak magnetic fields in the ISM, can inhibit the generation of 
Rayleigh-Taylor instabilities and/or help accelerate th e cool shell 
gas even further away through magnetic pressure (e.g. iFuiita et al] 
[2009,) . These effects influence the cold dense gas of bubbles, while 



the hotter, interior material is shown to escape to the hot halo in all 
of the simulations. Taking into account these results, we restrict the 
range of values of /r to /r ~ 1.1 — 2, implying that a significant 
fraction of the swept-up mass in bubbles stays in the ISM. The hot 
gas contained in the interior of bubbles is assumed to fully escape 
into the hot halo. In our standard model, we adopt /r = 1.5. In 
§ 14.3.21 we show how the mass outflow rate varies when /r takes 
the lowest and highest values in the range above. 

Fig. [2] shows a schematic of the evolution of bubbles in the 
ISM. We summarise all the parameters needed to characterise 
GMCs and the ISM of galaxies in Table|T] We give there the refer- 
ence value used for our standard SNe feedback model but also give 
the ranges motivated by observations and theory, which we also test 



3 INCORPORATING DYNAMICAL SUPERNOVA 
FEEDBACK INTO A GALAXY FORMATION 
SIMULATION 

One of the aims of this paper is to study how the outflow rate 
depends on galaxy properties in a galaxy population which has a 
representative set of star formation histories and which resembles 
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(1) Early stages of bubble evolution 




(2) Bubble at the point of breaking out of the ISM 




Figure 2. Geometry of the dynamical model for supemovae feedback. Top 
panel: the eai'ly stages of pressurised bubble growth due to SNe, where 
the bubble is fully embedded in the ISM. at a distance d from the galaxy 
centre, where the disk has a gas scale height of hg. The bubble radius and 
expansion velocity are and Vb, respectively. Bottom panel: Schematic 
showing the stage of bubble evolution just before breaking out from the 
ISM. At this stage the bubble has just exceeded the gas scale height. 



observed galaxy properties. We achieve this by incorporating the 
full dynamical model described in §|2]into the semi-analytic galaxy 
formation model GALFORM, which is set in the A cold dark matter 
framework. 

In §[3TT1 we briefly describe the GALFORM model and in ^ 13.21 
we give details on how we modify the model to include the dynam- 
ical model of SNe presented in ij|2land |2.2. 11 



3.1 The GALFORMmodel 

The GALFORM model takes into account the main physical 

S irocesses that shape the formation and evolution of galaxies 
Cole et alJl200Ql These are: (i) the collapse and merging of DM 
halos, (ii) the shock-heating and radiative cooling of gas inside 
DM halos, leading to the formation of galactic disks, (iii) quies- 
cent star formation (SF) in galaxy disks, (iv) feedback from super- 
novae (SNe), from AGN and from photo-ionization of the IGM, 
(v) chemical enrichment of stars and gas, and (vi) galaxy merg- 
ers driven by dynamical friction within common DM halos, which 
can trigger bursts of SF and lead to the formation of spheroids (for 
a review of these ingredients see Baugh 2006 and Benson 2010). 
Galaxy luminosities are computed from the predicted star forma- 
tion and chemical enrichment histories using a stellar population 
synthesis model. Dust extinction at different wavelengths is cal- 
culated self-consistently from the gas and metal contents of each 
galaxy and the predicted scale lengths of th e disk and bulge c om- 
ponents using a radiative tra nsfer model (see lLacev et al.l20T 3 and 
iGonzalez-Perez et al.ll2012h . 

GALFORM uses the formation histori es of DM halos as a start- 
ing point to model galaxy formation (see'Cole et al. 2000). In this 
paper we use hal o merger trees extrac ted from the Millennium N- 
body simulation jSpringel et al.ll2005l) , which assumes the follow- 
ing cosmological parameters: Qm = ^om + fibaryons ~ 0.25 



(with a baryon fraction of 0.18), Ha = 0.75, as = 0.9 and 
h — 0.73. The resolution of the A'^-body simulation corresponds 
to a minimum hal o mass of 1.72 x 10^"h^^ AIq, which in the 
iLagos et al.l ( 1201 2h model corresponds to a stellar mass limit of 
7 X 10^ Mq. This is sufficient to reso lve the halos that c ontain 
most of the H2 in the universe at z < 8 dLagos et alj|201 la ). The 
construction of the m erger trees used by GALFORM is described in 
lMersonetaLl l l2013h . 

In this paper we focuses on the Lagos et al. (2012; hereafter 
Lagosl2) model, which includes a two-phase description of the 
ISM, i.e. composed of the atomic and molecular contents of galax - 
ies, and adopt the empirical SF law of iBlitz & Rosolowskvl ( l2006h . 
The physical treatment of the ISM in the Lagos et al. model is a key 
feature affecting the predicted outflow rate of galaxies, as we show 
in § |4l which justifies our choice of exploring the full dynamical 
model of SNe in this model. 

The iBlitz & Rosolowskvl ( I2OO6I) empirical SF law has the 

form 



SsFR = i^SF fmol Sg 



(33) 



where Esfr and Egas are the surface densities of the SFR and 
the total cold gas mass, respectively, vsf is the inverse of the 
SF timescale for the molecular gas, usf = ^sf^, and fmoi ~ 
Enioi/Egas is the molecular to total gas mass surface density ra- 
tio. The molecular and total gas contents include the contribution 
from helium, while the HI and H2 masses only include hydrogen 
(helium accounts for 26% of the overall cold gas mass). The inte- 
gral of Esfr over the disk corresponds to the instantaneous SFR, 
tp. The ratio fmoi is ass umed to depend on the int ernal hydrostatic 
pressure of the disk as (iBlitz & Rosolowskvlliooel) 



— fmol/(fmol ^ 1) — 



Poxt 

"p7 



(34) 



For a description of how we calculate Poxt see Appendix IB II 
The parameter values we use for i^sf, Po and ap are the best 
fits to observations of nearby spiral and dwarf galaxies, vsf ~ 
0.5 Gyr-\ op = 0-92 and log(Po/kBfc r n~-^K]) = 4.54 
(iBlitz & Rosolowskvl |2003 : iLerov et aklboOSl : iBigiel et alj|201ll: 
iRahman et aDl2012h . 

For SBs the situation is less clear. Observational uncertain- 
ties, such as the conversion factor between CO and H2 in SBs, and 
the intrinsic compactness of star-forming regions, have n ot allowed 



a clear characterisation of the SF law in this cas e (e.g. iKennicutt 



19981: iGenzel et al.ll201(]|:ICombes et al1l201 ll: see lBallantvne et al 



20131 for an analysis of how such uncertainties can bias the in- 



ferred SF law). Theoretically, it has been suggested that the SF 
law in SBs is different from that in normal star-forming galaxies 
jPelupessv & Papadopouloil2009l) . The ISM of SBs is predicted to 
always be dominated by H2 independently of the exact gas pres- 
sure. For these reasons we choose to apply Eq. |33]only during qui- 
escent SF (i.e. SF fuelled by the accretion of cooled gas onto galac- 
tic disks) and retain the original SF prescription for SBs, which are 
driven either by galaxy mergers or disk instabilities (see lCole et al.l 
I2OOOI and Lll for details). In the SBs, the SF timescale is taken 
to be proportional to the bulge dynamical timescale above a min- 
imum floor value (which is a model parameter) and involves the 
who le ISM gas content in the SB, giving SF R = Mgas/TSF,SB 
(see lGranato et al.l2000l and lLacev et ai]|2008l for details), with 

TSF.SB = max(rmin,fdynT"dyn)- (35) 

Here we adopt T min = 100 Myr and fdyn ~ 50 following 
iLagos et al.l ( l2012l) . 
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Throughout the paper we will refer to galaxies as 'starburst 
galaxies' if their total SFR is dominated by the starburst mode, 
SFRstarbuist > SFRquicsccnt, while the remainder of the model 
galaxies will be referred to as 'quiescent galaxies'. 



ish the atomic/molecular gas contents and also modify the H2/HI 
ratio, as the gas and stellar surface densities change. 

The set of equations describing the flow of mass and metals 
between the different phases are 



3.2 Predicting the star formation history of galaxies 

The GALFORM model includes two gas phases in the ISM of galax- 
ies, an atomic and a molecular phase, which coiTespond to the 
warm and cold phases, respectively. By including dynamical mod- 
elling of SNe feedback, we introduce a new phase into the ISM of 
galaxies corresponding to the interiors of expanding bubbles (see 
§|2}. 

The equations of SF need to be modified accordingly to in- 
clude the contribution from the mass and metals in bubbles. The 
chemical enrichment is also assumed to proceed through the ex- 
pansion of SNe inflated bubbles: stellar winds and SNe feedback 
shock the surrounding medium and inflate bubbles through thermal 
energy, so the new metals produced by recently made intermediate 
and high mass stars will be contained in the interiors of bubbles. 
In the case of low mass stars, recycling of mass and newly syn- 
thesised metals feed the ISM directly. In the case of confinement, 
metals contained in the thin, dense shell of swept-up gas and the in- 
terior of bubbles are mixed instantaneously with the cold and warm 
ISM. Note that we do not apply any delay to the mixing of metals 
given that the cooling time for the hotter phases is typically small 
(tcooi = 5 X 10" - 10^ yr). 

The five mass components of the system are: the stellar mass 
of the disk. A/*, the total gas mass in the ISM (molecular plus 
atomic), MgjsM, the mass in bubbles (interior plus shell) in the 
ISM, Mb, ISM, the mass of the hot gaseous halo of the galaxy. A/hot, 
and the mass escaping the galaxy disk through bubbles, Meject- The 
latter represents all gas that has not yet mixed with the hot halo gas; 
i.e. that is thermally/kinematically decoupled from the hot halo gas. 
The underlying assumption is that all gas ejected from the disk ends 
up in reheated gas reservoir The reincorporation time, Trein, of the 
ejected component into the halo is always larger than the timestep 
over which we perform the integration. We therefore calculate the 
rate of reincorporation of gas into the hot halo component only with 
the ejected mass available at the beginning of the timestep, Afcjcct. 
We remind the reader that in this paper we we use the standard ap- 
proach of GALFORM to calculate r^d-n- This consists of parametris- 
ing Trcin as depending linearly on the dynamical timescale of the 
halo regulated by an efficiency, which is a free parameter of the 
model, Trcin = Tdyn/archcat (we retain the value of Qrchcat = 1.2 
used in Lagos 1 2). In paper II we introduce a physical modelling of 
the reincorporated gas and the timescale for this process. 

Fig. [3] depicts the exchange of mass and metals between the 
different components of galaxies: the hot halo, ISM, stars and bub- 
bles e xpanding in the ISM. As in the original model of ICole et al I 
( I200Q) . we assume that during SF, the inflow rate from the hot halo, 
Afcooi, is constant, implicitly assuming that SNe heating plays no 
role in the inflow rate until the ejected mass and metals are incor- 
porated into the hot halo after timescale Trcin. The gas mass in the 
ISM is affected by Afcooi, the rate at which mass is recycled from 
evolved stars (assumed to go straight to the ISM), the rate at which 
bubbles sweep up mass from the ISM, Afsw.iSM, and the rate of 
bubble confinement, A/conf,iSM and break-out, Afbo.iSM (the cal- 
culation of each of these are described in detail in Appendix [Ct. At 
each substep in the numerical solution scheme, we update the val- 
ues of each of the mass variables. It is therefore possible to replen- 



Mass exchange : 
Af* = (1 - i?ES - i?SN)V', (36) 
Afg,isM = Mcool + (-Res — l)'/' — A/sw.ism + Afconf.iSM 

+(1 - /bo)Afbo,iSM, (37) 

Afb,ISM = -RsN'i/' + A^sw.ISM — AfconfJSM " A^boJSM (38) 

Afcjcct = /boAfbojSM - (39) 
A^hot = -Afcooi + ^^. (40) 

Trcin 

Metallicity exchange : 

Aff = (l-i?ES-i?SN)^gV', (41) 
-^^ISM = Afcool^liot + (pes + i?ES^g)V' — Mg^jsM 

+ Mc^onf,ISM + (1 - /bo)A>4'^o,iSM, (42) 

,ISM 

(43) 

Mcfcct = /boAfb^o,iSM —, (44) 

"^rcin 

MLt = -Mcooi^hot + (45) 

Trcin 

The recycled mass from newly formed stars is specified sep- 
arately for SNe, Rsn, and intermediate and low mass stars, J?es 
(namely, evolved stars). We calculate the recycled fractions of each 
stellar mass range following Eq. IA2I SNe are considered to be 
all stars with m > 8A/q, and less massive stars in the range 
1 < m/MQ < 8 are considered as evolved stars (intermediate and 
low mass stars). Stars less massive than IMq have lifetimes larger 
than the age of the Universe and therefore do not recycle mass into 
the ISM. The yield is also defined separately for SNe and evolved 
stars in order to inject the metals from SNe into the bubbles, whilst 
metals from evolved stars go directly into the ISM. We adopt the in- 
stantaneous mixing approximations for the metals in the ISM. This 
implies that the metallicities of the molecular and atomic phases in 
the ISM are equivalent and equal to Zg = A/^jigi^/Airg,disk. The 
metallicity of the hot gas in the halo is Zhot = Mj^t/Afhot- 

The system of SF Eqs. 136143 1 applies for quiescent SF and 
SBs. In the latter case Afcooi = 0. During a SB, we assume that all 
bubbles expanding in galaxy disks are destroyed, as well as bubbles 
expanding in the satellite galaxy in the case of a galaxy merger. The 
new generation of stars made in the SB creates a new generation of 
inflated bubbles expanding over the bulge. 



4 PHYSICAL CHARACTERISATION OF BUBBLES IN 
THE ISM 

In this section we explore the physical properties of bubbles and the 
main drivers of their evolution in the ISM of galaxies. In ij |4. 11 we 
focus on individual examples of bubbles in ad-hoc galaxies. We ex- 
plore how the bubble mass depends on different global galaxy prop- 
erties, such as the gas fraction, gas metallicity and scaleheight, and 
local properties, such as gas density and surface density. In § 14.21 
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Figure 3. Schematic of the flow of mass and metals in the dynamical model of SNe feedback. The scheme shows the exchange of mass and metals (solid lines) 
between the hot halo, stars and the thi'ee gas phases in the ISM, and the partition of the ISM gas into the atomic and molecular gas components (dashed lines), 
corresponding to Egs. l36l4"3] in the text. Note that the same scenario would apply to SBs without the inflow of cooled gas from the hot halo. 
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Figure 4. Surface density of molecular and atomic gas as a function of the 
distance from the galactic centre in units of the half-mass radius of the three 
example galaxies listed in Table |2] Line styles and colours show different 
components of the gas content in the different galaxies as labelled. 



ij 14.31 and 14.41 we focus on the outflow properties of GALFORM 
galaxies when the full dynamical model for SNe feedback is in- 
cluded (see ij |2.2.2t . Comparisons with observations and previous 
theoretical work are presented and discussed in 14.41 



4.1 Properties of individual bubbles 

We study the dependence of the mass in a single bubble (interior 
plus shell) on the properties of the diffuse medium with the aim of 
determining which local properties are the more relevant in setting 
the mass of bubbles at the point of break-out or confinement (i.e. 
their maximum mass). 



Table 2. Properties of the three example galaxies used to study the effect 
of the different physical parameters on the evolution of bubbles in the ISM. 
We list the 10 properties we need to characterise the radial profiles of the 
stellar, gaseous and DM components, disk and bulge half-mass radii, 
and r^, stellar mass in the disk and the bulge, ^ and M^ ^,, cold gas 
mass, Mgas.iSM. gas metallicity, Zg, halo viiial mass, M^aio. radius, Tvir, 
and halo concentration, c. We also fix the distance to the galaxy centre at 
which the example bubble is located, d. The properties listed define the local 
properties of the ISM (see Appendix |B). For those parameters which we 
vary, we give the range chosen to study their effect on the bubble expansion, 
and in the Line below this we give the reference value. 



Model 


Dwarf 


Spiral 


Giant 


Varying pai'ameters 


Afgas,ISM/M0 


107.lo9-5 


108-10" 


lO^-lO^^ 


ref . value 
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1 X 10" 
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10^-1012 


ref. value 
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Figure 5. Top panel: Bubble mass at the point of break-out or confinement, maximum m^, as function of the distance to the galaxy centre in units of the 
half-mass radius, d/rrM (left panel), the gas metallicity in units of the solar metallicity, Z^/Zq (middle-left panel), the ISM mass (molecular plus atomic), 
Afg isM (middle-right panel) and the stellar mass, A/gtellar (right panel). The segments of the curves shown with solid lines correspond to those regions of 
the planes where bubbles end up breaking out from the galactic disk. Those segments shown with dashed lines correspond to regions where bubbles end up 
confined in the ISM of the galaxy. Bottom panel: Bubble mass at the point of break-out as a function of the local properties (i) atomic gas density, (ii) total 
(molecular plus atomic) gas surface density, (iii) surface density of total gas plus stars, (iv) gas scaleheight, (v) gas fraction and (vi) the ratio between the 
interior and the swept-up mass of bubbles (the interior mass corresponds to the fraction of the total mass injected by SNe that has not yet cooled down or hit 
the shell). Individual realizations for each galaxy are shown as points in the colours labelled. 



In order to fully characterise a single bubble in the ISM of a 
galaxy, we need to choose values for the galaxy properties which 
are required in the dynamical SN feedback model, namely the gas 
and stellar mass in the disk and the bulge, the half-mass radii of 
both stellar components, the halo virial mass, radius and concentra- 
tion, the gas metallicity and the location of the bubble in the galaxy 
disk. We focus on three example galaxies with properties within a 
representative range which are listed in Table |2] 

To calculate the expansion of a single bubble in the ISM of 
these galaxies, we use the standard set of parameters in Table [T] to 
describe GMCs and the ISM. In Fig. [4] we show the radial profiles 
of the atomic and molecular gas for t he three galaxies of Table [2] 
We construct these profiles using the iBlitz & Rosolowskvl ( l2006l) 
relation (Eg. |34t. The three galaxies plotted in Fig. |4] show cen- 
tral regions dominated by molecular gas, and atomic gas surface 
densities which saturate at ~ 10 Mq pc~^, above which the gas is 
mainly molecular. 

In order to study the dependence of the maximum mass of 
bubbles on galaxy properties, we vary the mass of gas and stars, 
the gas metallicity and the distance of the bubbles from the galaxy 
centre for the three galaxies in Table |2] These parameters are ex- 
pected to have an effect on the expansion of bubbles by varying 
the gas density, scale height, cooling timescale, gravitational field, 
etc. The strategy is to vary one property at a time leaving the other 
ones unchanged, to see how the predictions change. We evolve bub- 
bles until they become confined or break out from the galaxy disk. 
When we fix d, we arbitrarily choose d = 0.5 rso for illustration. 
This value of d typically coiTesponds to a region where bubbles 
break out. The 4 experiments (i.e. changing d, Z^, A/gas, ism and 
Mi,,d) are performed for each of the galaxies of Table |2] and the 



results are shown in the top panel of Fig.|5] The maximum mass of 
a single bubble shown in Fig.|5]corresponds to the mass at the point 
of break-out or confinement. 

In the central regions of galaxies, bubbles break-out from the 
galaxy disk, while in the outskirts bubbles tend to be confined. In 
the case of the 'dwarf galaxy, the break-out region is restricted to 
d < O.Srso, while in the case of the 'spiral' and 'giant' galaxies, 
the region of break-out extends out to d > rso . In the break-out 
regions, there is a strong relation between the bubble mass and the 
distance from the galactic centre. This is driven by an underlying 
relation between mb and the gas scaleheight or gas surface density. 

Variations in the gas metallicity have very little effect on the 
resulting bubble mass. When the gas surface density is high, the 
metallicity plays only a minor role because the cooling time is al- 
ready very short and bubbles become radiative very quickly. In the 
case of low gas surface densities, the cooling time becomes long 
even for high metallicities, which preserves the energy of the bub- 
bles. In the case that metallicity does have an effect on the bubble 
mass, the differences found are always less than a factor of ~ 2. 

Strong variations in the maximum mass of the bubble are ob- 
tained when varying Mgas,iSM- In the regime of break-out from 
the galaxy disk, the bubble mass quickly decreases when increas- 
ing A/gas, ISM- As Afgas.iSM Increascs, the surface density of gas 
also increases. This reduces the gas scaleheight, which reduces the 
bubble mass. The reason for this is that the radius the bubble needs 
to reach to escape the galaxy decreases, and therefore also the to- 
tal mass that it is able to sweep-up also decreases, as this is pro- 
portional to the bubble volume. The higher Afgas,iSM results in an 
overall decrease of the bubble mass by a factor of 100 — 500. 

Variations in stellar mass have a non-negligible effect on the 
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bubble mass, particularly at the massive end of the range tested 
(see second row of Table |2ll. There is a trend of decreasing bub- 
ble mass with increasing stellar mass in the region of break- out. 
This happens due to the increasing gravitational field driven by the 
higher stellar surface densities, which decreases the gas scaleheight 
of the disk and the radius the bubble needs to reach to break-out. 
The bubble mass obtained when increasing the stellar content of 
galaxies can be lower by up to a factor of 3. The effect of the more 
efficient deceleration of bubbles due to the larger gravitational field 
when the stellar mass increases is secondary to the effect of the 
stellar surface density on the gas scaleheight, and represents only 
~ 0.1 — 5% of the total effect observed when increasing Al^^d- 

The distance to the galactic centre and the gas content of the 
galaxies shown in Fig. |5] drive the strongest variations in bubble 
mass. This is due to the dependence of mb on the gas density 
(atomic plus molecular) and the gas scaleheight, which is shown 
in the bottom-right panel of Fig. |5] We include only those exam- 
ples in which the bubble breaks out from the galaxy disk. Bubble 
masses in the cases tested here are always dominated by the swept- 
up mass (see bottom-right panel of Fig. |5}- However, there is an 
increasing contribution from mint to mb for decreasing mb. We 
give physical insight into the relations between mb, /ig and Egas in 
the next subsection. 

In the case of the gas fraction, we find that there is a complex 
dependence of mb on /gas- The gas fraction acts to modify the 
normalisation of the relation between the total outflow rate and /ig 
and the power-law index of the relation between the total outflow 
rate and E^as. 



70Mq pc ^, the power-law relations mb oc and mb cx 

Sjas , and in the lower density regime, we find mb oc h°g'' and 
mb oc Egas**- These power-law relations are approximate as the 
exact value of the power-law index changes slightly from case to 
case. From this analytic derivation of the scaling relations it is fair 
to say that the transition from the atomic- to molecule-dominated 
media has a large impact on the mass of a bubble at the point of 
break-out. 

If we assume a steady state (i.e. the SFR is constant), we can 
write the outflow rate per annulus as a function of each individual 
bubble mass as. 



/bo rUh Mmol 
Tlifo.GMC MqmC 



(50) 



Considering tp = i^sf Afniob we can directly write /3 per annulus 
in terms of a single bubble mass 



a Afojcct 
P = 



/b 



i^SF Tlife,GMC MgMC 



(51) 



There is a direct relation between /3 and mb in the case of a steady 
state. We therefore expect to see a similar transition in the relation 
between the outflow rate and the gas surface density to the one ob- 
tained for mb: from a steeper relation in galaxies with molecule- 
dominated ISM to a shallower relation in galaxies with atomic- 
dominated ISM. From Eqs.|48]and [5T]we also see how each of the 
parameters describing the ISM and GMCs affect individual bubble 
masses and the global outflow rate. 



4.1.1 Analytic derivation of the scaling relations of single 
bubbles 

At the point of break-out, the volume of the gas disk occupied by 
a single bubble is V = 27: hg{f^ — 1/3). In the regime where 
minj <C msw, which is a representative limit for most bubbles (see 
the bottom-right panel of Fig.|5}, and neglecting temporal changes 
in the gas density of the diffuse medium during the evolution of 
bubbles in the ISM, one can write the bubble mass as 



mb = Pd V = (1 - /mol)7r(/,. - 1/3) Egas ftg 



(46) 



In order to find an expression for mb in terms of Egas and /ig alone, 
we need to express /moi as a function of the same variables. 

We can write f^oi in terms of the gas (atomic and molecular) 
density 

1 1 



l + (Pe./P0)- l+(^aI/Po 

By introducing the expression for /moi into Eq.|46l we find that 



(47) 



mb 



7!" (/r - I) Sg h, 



1\ ( 2-Po 
t1 — a 



■ al/Po \ < 1 
■'^I/J^o)»l (48) 



If we now apply the limit Eg 3> (crg/(Tt)E*, where gas dominates 
over stars in the gravity acting on the gas layer, we find 



mb oc 



ftg a Eg 
/ig+^"P oc E, 



-(l + 2ap) 



.fmol < 1 
.fmol ~ 1 



(49) 



These expressions describe the relations shown in the bottom panel 
of Fig |5] where we obtain, in the high-density regime, Egas ^ 



4.2 Radial profile of the mass loading factor and outflow 
velocity 

In order to physically characterise the outflow rate in a galaxy pop- 
ulation which resembles the observed one, we use the GALFORM 
semi-analytic model, into which we incorporate the dynamical 
feedback described in § [2] The key difference with the analysis of 
§|4T|is that here we explore the whole galaxy population and the 
outflow rate with the aim of characterising: (i) a preferred radius 
from which most of the material escapes and the outflow velocity, 
and (ii) the scaling relations between the mass loading factor, /?, 
and local properties of the disk, computed in an annulus which is at 
a distance d from the galactic centre. The galaxies used in the anal- 
ysis in this section are selected so that they are close to the break of 
the stellar mass function at low-redshift, > 10^" Mq , and 
have 2 < 0.1. This selecti on makes the galaxy properties compa- 
rable to those simulated bv lCreasev et al.l fcoiSh . 

In order to gain insight into (i), we show in the top panel 
of Fig. [6] the outflow rate in each radial annulus in units of the 
global SFR as function of the distance from the galactic centre. We 
distinguish between galaxies with different gas fractions, /gas = 
Afg,isM/(AfgjsM + M*). There is a tendency for gas-rich galax- 
ies to have most of the mass breaking-out from the disk at d ~ rso, 
while in gas-poor galaxies most of the mass escapes from close 
to the galactic centre. We calculate the radius inside which half of 
the global outflow mass escapes, MovLt{d < rout) = Mcjoct/2, 
where A/cjcct is the global outflow rate. Galaxies in Fig. |6] with 
/gas > 0.8 have rout = 0.8 rso and those with /gas < 0.1 have 
rout = 0.4 rso. This is consistent with the picture presented in 
§ 14.11 where the gas-poor dwarf galaxy has a more centrally con- 
centrated outflow than galaxies that are gas rich. Given that high- 
redshift galaxies are pred icted to be more gas-rich than their lo w- 
redshift counterparts (see lGeach et al.EoilHLagos et al.ll201Iah . it 
is therefore expected from this model that high-redshift galaxies 
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Figure 6. Top panel: The outflow rate contributed by each annulus in units 
of the global SFR, as a function of the distance from the galactic centre 
in units of the half-mass radius, d/r^g, for the dynamical model with the 
standard choice of parameters (see Table [T] and for galaxies at z < 0.1 
and with > 10^" /i~^Mq. For quiescent SF we use rso of the disk, 
and for starbursts, rso of the bulge. Sohd lines and en'orbars represent the 
median and 10 to 90% range of the distributions. The predictions are plotted 
for galaxies with different gas fractions, as labelled. Bottom panel: As in the 
top panel, but here we show the outflow velocity of the gas at the point of 
break-out as a function of distance from the galactic centre in units of the 
half-mass radius. 




12 3 4 

log(Pgas+2:-]/MePC^) 

Figure 7. Top panel: The ratio of the outflow rate to SFR per annulus as a 
function of the surface density of gas plus stars for galaxies at z < 0.1 and 
with Mt > 10^*^ h~^MQ and for different gas fractions, as labelled, in 
the model with the standard set of parameters (see Table[T]. Solid lines and 
errorbars coiTespond to the median and 10 and 90 percentile s of the distri- 
butions. The shaded region corresponds to the predictions of lCreasev et alj 
(2013), and is plotted over the range of surface density of gas plus stars 
probed by the simulations. Bottom panel: As in the top panel but for the 
outflow velocity per annulus as a function of (Eg + S«). 



exhibit more extended outflows relative to the size of the galaxy 
than low-redshift galaxies. 

In the bottom panel of Fig. [6] we show the mass-weighted 
velocity of the gas escaping the galaxy disk as a function of 
the distance from the galactic centre, d, for galaxies with differ- 
ent gas fractions. There is a trend of increasing outflow velocity 
with increasing /gas. Gas rich galaxies typically have a molecule- 
dominated ISM. In these galaxies the density of atomic, diffuse 
gas is lower, resulting in a more inefficient deceleration of bubbles. 
The predicted values of the outflow velocity are comparable with 
the observed values. We directly compare with observations of the 
outflow velocity in i; 14.41 

Concerning the scaling relations of the outflow (listed as (ii) 
above), we calculate the ratio between the mass outflow rate and 
the SFR in each annulus, /3annuius, and investigate its dependence 
on the local properties of the disk, as estimated at the mean radius 
of each annulus. The top panel of Fig.[7]shows the relation between 
/3annuius and (Eg + E*), evaluated at rannuius, for galaxies with 
different gas fractions. There is a tight correlation between the two 



quantities, with only a modest dependence on other galaxy proper- 
ties, such as the gas fraction. This is expected from t he correlation 
betwe en mb and (Eg + E*) (ij 14.1b . The results of ICreasev et al.l 
( l2013h (see §[T]for details) are also shown in Fig.|7]by the shaded 
region, plotted over the range of surface densities probed by their 
simulations. Our predicted relation is similar to what Creasey et al. 
found using a completely different approach (see §[T). 
The best fit to the relation in Fig.|7]is 



annulus 



+ E„ 



69M0PC- 



(52) 



The bottom panel of Fig. [7] shows the outflow velocity, 
^outflow, as a function of (Eg + E^,), evaluated at rannuius. There is 
a trend of increasing ^outflow for increasing (Eg + E^,). Our predic- 
tions for ^outflow also overlap with those of Creasey et al., although 
we find that outflow velocities > 1000 kms^^ are statistically un- 
likely. These velocities can occur for starbursts in our model (see 
§|4XTJ. Note that for a given (Eg + E^,) there is a trend of 13 de- 
creasing with and ^outflow increasing with increasing gas fraction. 
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This prediction is also in agreement with the findings of Creasey et 
al.. 

Note that changes in the SNe feedback model parameters, 
which are summarised in Table [T] produce similar deviations to 
those found for the galaxy-wide /3 and mass-weighted ^outflow in 
§ 14.3.21 We find that the surface density normalisation and power- 
law index in Eq. [52] increase with increasing redshift, in a similar 
way that the global /3 does (Fig. 114b. Therefore, the similarity be- 
tween our predictions and those of Creasey et al. is confined to our 
low-redshift galaxy sample. Note that the results of Fig. [7] for a 
fixed gas fraction do not depend on stellar mass or redshift, but the 
global normalisation and power-law index of Eq. |52]do due to the 
predominance of gas poor galaxies at low redshift and of gas-rich 
galaxies at high redshift. 



4.3 Statistical properties of tlie outflow rate and velocity 

In this section, we attempt to answer three questions: What is the 
effect of the multiphase treatment of the ISM on /3? What is the 
overall effect of varying the physical parameters of the ISM and 
GMCs on the outflow rate? Is the outflow rate dominated by adia- 
batic or radiative bubbles? 

Here we analyse galaxies from GALFORM, after the full dy- 
namical model of SNe feedback is included in the calculation. 
At each redshift we focus on galaxies with A/* > 10* h~^MQ, 
to be safely above the resolution limit of the Millennium simula- 
tion (§ I3.lt . We consider the total mass loading rate of the out- 
flow, P, which we define as /3 = iV/cjcct/V'> where Afejoct cor- 
responds to the total mass breaking out from the ISM (given by 
/boA/bo.isM in Eqs |36|43t and i/; is the instantaneous SFR. In 
ii l4.3.4l we analyse the metal loading of the wind, which we define 



as 13^ 



Af4,t/Zgt/). This /3 differs from the /3annuius of §|42] 



in two respects; the former is integrated over the galaxy and over 
longer timesteps. 

In ij |4.3.1ll4.3.3l and l4.3.2l we show the total mass loading /3 as 
a function of the gas scaleheight at the half-mass radii of galaxies, 
hg. This can be understood from the strong dependence of mb on 
hg and the small dispersion in this relation (see ii l4.Ib . In ij |4.3.4l we 
show how and where differs from /? and the reasons for such 
differences. 



4.3.1 Testing the ejfect of the multiphase medium and gravity on 
the outflow properties 

The top panel of Fig. [8] shows the correlation between (3 and hg at 
the half-mass radius obtained with and without considering gravity 
from stars and DM in Eqs. |I|3l|18|20| and l24l26] and using the stan- 
dard set of parameters to describe GMCs and the ISM of galaxies 
(see Table [TJ. We find that j3 is only slightly affected when grav- 
ity is not included. This agrees with what we find for individual 
bubbles, in which gravity has an effect of at most 5% on the final 
bubble mass. 

The effect of including the H2/HI ratio calculated from the 
Blitz & Rosolowsky pressure law in the modelling of the ISM is 
much larger than the direct gravitational effect, as the dotted line in 
Fig.[8]shows. The omission of self-consistent multiphase modelling 
is represented by the results obtained with a fixed H2/HI= 0.37 
ratio, which is the value used in p revious w ork to esti mate HI 
from the total cold gas content (e.g. IPower et a l. 2010; Kim et al.l 
I20T Ij). With a fixed H2/HI ratio, the mass loading increases by fac- 
tors of up to 100 for galaxies with the smallest gas scaleheights 
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Figure 8. Top panel: The mass loading, /3, as a function of the gas scale- 
height at the half-mass radius for: quiescent (solid line), starburst (dashed 
line) and a subsample of massive galaxies. Mi, > 10^° h'^ Mq (long- 
dashed line), in the model with the standard set of parameters (Table [T). 
In the case of quiescent SF, hg is evaluated at rso of the disk, and for star- 
bursts, at rso of the bulge. We include in the plot all galaxies in GALFORM at 
z < 1 and with Ah > 10** /i~^Mq. We also show the effect of suppress- 
ing gravity on the expansion of bubbles (dot-dashed line), and of assuming 
a constant H2/HI ratio instead of that derived from the Blitz & Rosolowsky 
pressure law (dotted line). Solid lines and errorbars indicate the median and 
10 and 90% ranges of the predictions. For clarity, errorbars are shown only 
for selected cases. Bottom panel As in the top panel, but here we show the 
mass-weighted outflow velocity as a function of the gas scaleheight. 



(i.e. highest density regi mes). This is due to the anticorrelation 
between H2/HI and hg l lLagos etalj|2011ah . Galaxies with very 
high gas and/or stellar surface densities have smaller hg and larger 
H2/HI, driving a lower overall content of HI and therefore provid- 
ing less material for bubbles to sweep up, reducing the outflow 
mass. This effect is very large in more extreme cases, where the 
pressure law predicts little HI. This is also clear from the single 
bubble examples of § 14.11 in which the bubble mass is greatly re- 
duced in molecule-dominated media. This demonstrates the impor- 
tance of the ISM m odelling introduced in lLagos et al.l ( 1201 l b) and 
iLagos et"aL ( 201 lah . and also included in some other recent models 
(e.g. lFuetalJl2010h . 

In the top panel of Fig. [8] we show the relations for starburst 
and massive galaxies separately. This stresses the similarity be- 
tween the relations displayed by quiescent and starburst galaxies in 
the P-hg plane and the fact that massive galaxies follow the same 
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relation as the overall galaxy population, which is dominated in 
number by lower mass systems. This is because the mass loading j3 
is primarily determined by the gas scaleheight and the gas fraction, 
as we show later in i) l5.2l 

In the bottom panel of Fig.[8]we show the mass-weighted out- 
flow velocity as a function of the gas scaleheight. There is a trend 
of decreasing velocity for increasing h^,. Starburst galaxies exhibit 
a relation with a similar slope to that of quiescent galaxies but off- 
set by « O.Sdex to larger velocities. This is due to the different star 
formation laws assumed in the model for the starburst and quies- 
cent star formation modes (see I3.lt . For a fixed h^, a starburst 
galaxy generally has a larger SFR than its quiescent counterpart. 
This drives larger energy and momentum injection, resulting in 
larger outflow velocities. However, the effect is not necessarily lin- 
ear, as the bubble deceleration also depends on the expansion speed 
(see Eqs.[3l[T8landl24t. The effect of gravity in the outflow velocity 
is only minor, as is also the case for /?. The effect of including the 
Blitz & Rosolowsky pressure law in the modelling of the ISM on 
the outflow velocity is more significant, and its omission results in 



velocities that are larger by a factor of ~ 2 at small h^. In i; 14.41 we 
compare our predicted velocities with observations. 



4.3.2 Assessing the impact of ISM and GMC parameters on the 
outflow properties 

The top panel of Fig.|9]shows the predicted mass loading as a func- 
tion of the gas scaleheight when varying the parameters associated 
with the modelling of GMCs and the diffuse medium (see TableO. 
Changes in the GMC and diffuse medium model parameters drive 
different normalisations in the /J-Zig relation but have a weak im- 
pact on the shape of the relation. The variations between the models 
that produce the smallest and largest [3 values, which correspond to 
adopting /r — 1.1 and vqy ~ 0.3 Gyr~^, respectively, are at most 
a factor of ~ 10. It is reasonable to argue that a better understanding 
of the multi-phase nature of the ISM and the properties of GMCs is 
very important, even more so than including some of the physical 
mechanisms in the expansion of bubbles, such as gravity. This was 
also hinted at in Fig. [8] from the effect of adopting a multi-phase 
ISM description of the outflow rate. 

The effect of each of the parameters in Table [T] on /3 is sum- 
marised below. 

• Smaller values of /r result in smaller /3 values by a factor ~ 
3 — 5. This is expected from the role /r plays in determining the 
break-out radius of bubbles and therefore the bubble mass (Eo.l48b. 

• Adopting a smaller SF coefficient or a smaller GMC mass 
drives an increase in /3 due to the anticorrelation of /3 with both 
usF and Mgmc (Eq.lSU. The effect of increasing vsf or A/gmc 
is therefore a smaller /3. Adopting a longer lifetime for GMCs also 
decreases /? due to the anticorrelation between /3 and riitc,GMC- 

• A smaller hydrostatic pressure normalisation in the Blitz & 
Rosolowsky law (see §[3]( drives larger /3 but only in galaxies which 
have a molecule-dominated ISM, as it only affects this regime (see 
Eq.|48b. In these cases, the lower Po drives smaller individual bub- 
ble masses and therefore smaller /? (see Eq.lSIt. Similarly, the ef- 
fect of decreasing dd is to slightly decrease P, which is also ex- 
pected from the analysis of ij 14.1. 11 

The effect of varying the parameters above on the mass- 
weighted outflow velocities, ^outflow, is shown in the bottom panel 
of Fig. 121 Variations in ^outflow due to different ISM parameter 
choices are smaller than in the case of /3, with a difference between 
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Figure 9. Top panel: The predicted mass loading, /3, as a function of the 
gas scaleheight. In the case of quiescent SF, /ig is evaluated at rso of the 
disk, and for starbursts, at rso of the bulge. The predictions ai'e shown for 
different choices of the model parameters, as labelled. We include in the 
plot all galaxies in GALFORM at z < 1 with Af* > 10* /i~^AfQ. Lines 
and errorbars indicate the median and 10 and 90 percentile ranges of the 
relations. For clarity, the percentile range is shown only for one model as 
they are all similar. Solid lines are used for the model with the standard set 
of parameters and those predicting the lowest and the highest /3 for a given 
hg. Dashed lines aie used for the rest of the models (see Table [T]. Bottom 
panel As in the top panel, but here we show the mass-weighted outflow 
velocity as a function of the gas scaleheight. 



the minimum and maximum ^outflow of ~ O.Sdex. The models pre- 
dicting the highest and lowest /3 are not the same as those predicting 
the highest and lowest ^outflow This is because iioutflow is more af- 
fected by those parameters directly changing the energy injection 
into the ISM by SNe. Indeed, the parameter that is most important 
in setting ^outflow is the star formation coefficient, vsf. The more 
efficient the conversion from gas to stars, the higher the outflow 
velocity. This is consistent with what is shown for quiescent and 
starburst galaxies in Fig. [7] 



4.3.3 The physical regimes of the outflow 

Bubbles inflated by SNe feedback can escape the galaxy in any of 
the three evolutionary stages described in § |2] We now quantify 
where and when each of these stages dominates the outflow of ma- 
terial. 

Fig. [To] shows the mass loading, /3, as a function of the gas 
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Figure 10. The mass loading factor, f}, as a function of tiie gas scalelieight at 
the half-mass radius for galaxies with Af* > 10** h'^Mp) in three redshift 
ranges, as labelled in each panel. In the case of quiescent SF, /ig is evaluated 
at rso of the disk, and for starbursts, at rso of the bulge. The contribution 
to the total /3 (solid line) from bubbles escaping in the adiabatic, pressure- 
driven and momentum-driven snowplough phases are shown as dashed, dot- 
dashed and dotted lines, respectively. The ratio between the rate of mass 
confinement and the SFR, /3confi is shown as triple-dot-dashed line. Lines 
represent the medians and the eiTorbars, which are shown for claiity only 
for the total /3, represent the 10 to 90 percentile range. 



scaleheight, /ig, evaluated at the half-mass radius for the model 
with the standard set of parameters. We find that at high redshift, 
most of the outflow in galaxies is produced by bubbles escaping in 
the momentum-driven stage, while low-redshift galaxies with small 
gas scaleheights have mass outflow rates dominated by bubbles es- 
caping in the pressure-driven stage. High-redshift galaxies have a 
gas scaleheight set by the gas surface density with a negligible 
contribution from the surface density of stars. In the low-redshift 
regime, galaxies with small gas scaleheight have, by comparison, 
a more important contribution from the stellar component. In fact, 
the median gas fraction of the galaxy sample with /ig < 10 pc in 



the high- and low-redshift samples is 0.98 and 0.18, respectively. 
Galaxies which have the gas scaleheight set mainly by the stellar 
surface density, have bubbles where the cooling time for the inte- 
rior gas is large enough for bubbles to escape the disk in the pds 
stage. In the case of the larger gas scaleheight galaxy population, 
the scaleheight set mainly by the gas surface density, so no signifi- 
cant difference with redshift is obtained. 

When bubbles escape the ISM in the radiative phase (i.e. pds 
or mds), this implies that most of the outflow mass is in a cold, 
dense phase (i.e. molecular or neutral atomic gas) and that the in- 
terior mass of the bubbles is only a minor contributor This qual- 
itatively agrees with what is observed in local galaxies (e.g. Tsai 
et al. 2012a,b). A quantitative comparison will be presented in a 
forthcoming paper (Lagos, Baugh & Lacey, in prep.). 

The adiabatic phase only rarely dominates the outflow rate, 
since the transition from the ad to the pds stage takes place early 
on in the evolution of bubbles. This transition almost always takes 
place on a timescale of ~ 10'^ — 10^ years. Full confinement due 
to deceleration of bubbles rarely takes place (i.e. the case in which 
no bubbles break-out from the galaxy disk), and happens mainly 
in places where the scaleheight is large and the bubble has time to 
decelerate to the velocity dispersion of the diffuse gas (i.e. at low 
gas densities). Most of the gas which remains in the ISM therefore 
corresponds to gas expanding in the direction close to the plane 
(i.e. the fraction (1 — /bo) in Egs. 136143] ) rather than to bubbles 
which are fully confined in the ISM. The tendency we find for bub- 
bles t o break-out in the radiative phase contrasts with what lMonacol 
l;2004b) found, whose model predicts that most bubbles escape dur- 
ing the adiabatic phase. This difference may be due to the assump- 
tions Monaco makes that bubbles expand against the hot phase. In 
our model, bubbles expand against the warm phase, whose density 
is typically higher than the hot phase, which results in larger cool- 
ing rates. We find that our approach gives answers more similar to 
fully hydrodynamical simulations in the range where they overlap 
(see §|tl. 



4.3.4 Outflow rates of mass and metals 

We have analysed the physics behind the dependence of /? on 
galaxy properties and gave analytic derivations for such relations. 
However, a key part of the impact of outflows on galaxy evolution 
is the fate of the metals carried away by bubbles. In the model, we 
assume that the metals which flow out from the galaxy accumulate 
in the ejected mass component, which is later reincorporated into 
the hot halo gas (see Eas. l36l43] l. The amount of metals outflowing 
from the galaxy therefore has a direct impact on the cooling rate 
of the hot halo gas and hence on subsequent gas accretion and star 
formation in the galaxy. 

Here, we analyse the loading factor of metals defined as — 
M^icct/{Zgip) (see Eq. |44ll. The top panel of Fig. [TT] shows the 
metal loading factor as a function of the mass loading factor for 
galaxies at different redshifts. Galaxies al z < 2 follow a relation 
which is close to — (3, but which shows a flattening at /3 < 0.5 
(i.e. in the small gas scaleheight regime). However, as the redshift 
increases, deviations become important and begin at increasingly 
larger /3. At z > 6 there is almost no correlation between and 
P, with ~ 30 independent of /?, albeit with a large dispersion. 
This behaviour is due to high-redshift galaxies having intrinsically 
lower metallicity gas from which stars form. In the low-metallicity 
regime, metals in bubbles coming from the swept-up gas are neg- 
ligible compared to those coming from SNe ejecta; in the limit of 
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Figure 11. Top panel: The metal loading factor, / 

Mo: 



/(ZgV>), as 
t/t/", for galaxies with 



a function of the mass loading factor, /3 

Mi, > 10* h~^MQ. Both quantities are integrated over the galaxy and 
in the same timesteps. Lines with errorbars represent the median and the 
10 to 90 percentile range, respectively, for galaxies at different redshifts, as 
labelled. The thick, straight line shows f}^ = fi. Bottom panel: normalised 
distribution of /3 for galaxies in the same redshift ranges as in the top panel. 



Zg < psN and AttR ZgPd«s < 
outflow rate due to a single bubble as 



we can write the metal 



''eject 



/bo PSN V'GMC = /bo PSN VSF MgmC ■ 



(53) 



The rate of metals flowing out from the galaxy in a given annu- 
lus is regulated by the number of GMCs in that annulus M, 



eject 



/boPSN yS¥ Mmol- 

regime 



We then calculate j3 per aimulus in this 



/bo PSN 



Zf, USY Mn 



(54) 



Because we assume instantaneous mixing in GALFORM, this /3 
is representative of the global metal loading factor In the limit of 
Zg <C PsN, shows no dependence on hg. However, the mass 
outflow rate has a strong dependence on Egas, regardless of the 
metallicity of the ISM. This results in very little correlation be- 
tween and /3 in this low-metallicity regime. 

If the ISM is already enriched with some metals, which coiTe- 
sponds to approximately Zgas > 0.05 — O.IZq, the density of the 
gas in the ISM also has an important effect on (3^ given that the 
term -iiTR^ ZgpdVe becomes comparable to or larger than the term 
rhj^-ij in the evolution of single bubbles (see Eq[3](. In this case, a 
correlation between and /3 arises. 

Although a non-linear relation between /3 and is predicted, 
we find that most galaxies in our simulation follow a relation which 
is close to — p. This can be seen from the distribution of /3 for 
different redshifts in the bottom panel of Fig.[TT] Quantitatively, at 
least 75% of galaxies at any redshift have /3 > 1 and at least 50% at 
z < 5 have /3 > 10. This puts at least half or more of the galaxies 
in the regime where ~ /3. Galaxies deviating this relation are 



the most metal-poor ones, which typically coiTespond to those with 
low stellar masses. As we show later in §|6l the inclusion of a metal 
loading factor with an independent parametrisation from the mass 
loading factor in GALFORM, has a small effect on the luminosity of 
galaxies. However, if we wish to analyse in detail the gas content 
of galaxies and the evolution of the mass-metallicity relation, we 



would need to allow for such variations in the /? 
included in the model. 



parametrisation 



4.4 Comparison with observations and non-cosmological 
hydrodynamical simulations 

We compare our predictions for the mass loadin g of th e wind , 
j3, with the values infe rred from observation s by Martini (Il999h . 
iMartin et al.l ( I20l3) and lNewman et al] ( l20I2h . who use absorption 
features in galaxy spectra (directly prob i ng wa rm gas in the cir- 
cumgalactic medium), and lBouche et alj ( 1201 2h . who use absorp- 
tion lines in the lines-of-sight to background quasars (prob i ng the 
outflow and inflow of gas). Martin (1999,) and .Bouche et al.l(l2012l) 
focus on L* galaxies at low redshift {z < 0.1), while Martin et al.l 
| |2012[) focus on galaxies at z « 1 and Newman et al. on galaxies at 
z « 2. In the case of Martin et al., we focus on the subset of galax- 
ies with multiwavelength data for which measurements of stellar 
mass and SFR were possible and those which have an outflow de- 
tected at the 3a level, or better. 

The top panel of Fig.ll2lshows /3 as a function of stellar mass 
for our standard model (see Table [T). Symbols show the median 
stellar mass and the jS inferred from observational samples. Error- 
bars on the mass axis represent the range of stellar masses in the 
observed galaxy populations and on the y-axis, the typical errors 
in the inferred /?. Our model predicts /3 values which are in good 
agreement with those inferred from observations. However, there 
are large uncertainties associated with the inference of outflow rates 
from observations, in addition to the statistical uncertainties aris- 
ing from the small number of objects sampled. The main uncer- 
tainties in the calculation of outflow rates from observations come 
from the conversion between the ion and hydrogen column densi- 
ties, which depends on the gas metallicity and ionization factor, and 
the assumed geometry in the case of classical absorption line spec- 
troscopy (e.g. IProchaska et ai]|20I ih . and the still uncertain nature 
of absorption by low-ionisation metal lines, in the case of absorp- 
tion line studies in the quasar sightlines. 

The bottom panel of Fig.[T2lshows the mass-weighted outflow 
velocity as a function of stellar m ass. We show the observational in- 
ferences of lBouche etai] ( l20I2h at 2 « 0.1 and lMartin et alj ( l2012h 
at 2 1. In the case of Bouche et al., Mgll absorption lines are 
used to infer an outflow velocity, while in the case of Martin et 
al.. Fell absorption lines are used. The predicted outflow veloci- 
ties agree with those inferred from the observations, although the 
physical interpretation of the inferred velocities and outflow rates 
is not straightforward, as the different gas phases of the outflow 
could have different velocities and mass loadings. In the case of 
the model, the plotted outflow velocities are calculated from the 
expansion velocities of bubbles at the point of break-out and are 
dominated by the phase that contributes the most to the outflow 
mass. We predict that in many cases this corresponds to a warm or 
cold phase (neutral or molecular). In the case of observations, the 
available data probe warm ionised gas and are corrected to account 
for the neutral component. Ideally, these data need to be comple- 
mented by deep observations at millimeter wavelengths to directly 
probe the part of the outflow that is in a cold phase. In a future 
paper we analyse more fully the outflow mass in different phases 
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Figure 12. Top panel: The mass loading, /3, as a function of stellar mass 
for galaxies in three different redshift ranges for the standard set of parame- 
ters (Table[T]. Solid lines and eiTorbars indicate the median and 10 and 90% 
ranges of the distributions. For clarity, en'orbars are shown only for the low- 
est and highest redshift ranges. The observationally inferred /3 are shown 
using symbols. These data are from J VIartin ( 1999) at 2; O.l. lMartin et alj 
i2012l) aX z 7^ 1 and lNewman et al.N2012l) at z ^ 2 using galaxy ab- 
sorption line spectroscopy and from lBouche et al.] ( l2012l) at z 0.1 using 
spectroscopic studies of lines-of-sight to background quasars. The errorbars 
in the mass axis for the observational points represent the range of stellar 
mass in each sample and in the i/-axis we show representative eiTors in the 
inferred /3. Bottom panel: The mass-weighted outflow velocity, tJoutflow. 
as a function of stellar mass for galaxies in three different redshift ranges 
for the standard set of parameters, as in the top panel. Solid lines and er- 
rorbars indicate the median and 10 and 90% ranges of the relations. For 
clarity, eiTorbars are shown for the lo west and highe s t reds hift ranges. The 
observationa l ly infe rred /3 values from lBouche et al] j2012l) at z 0.1 and 
lMartinetal]j2012l) at z ^ 1 are shown as symbols. 



and carry out a more detailed comparison with observations (La- 
gos, Baugh Sl Lacey, 2013b). 

We find that our model agrees better with observationally in- 
ferred outflow rates compared to previous work on SNe feedback 
and mass ejection from the ISM. For example, lEfstathioul fcOOOh 
implemented a physical model for galaxy evolution in which self- 
regulation was imposed: energy loss by cloud collisions is compen- 
sated by the energy input by SNe. Efstathiou predicted that galax- 
ies with Mateiiar ~ 5 X 10^" Mq havc a mass loading factor in 
winds from the ISM of /3 ~ 0.2, which is a factor of more than 10 



lower than the values inferred bv lMartinl(ll999h and lSouche et al.l 
| |2012|) . The assumptions in the modelling of Efstathiou are differ- 
ent from ours. An important difference is that we do not assume 
self-regulation in galaxies but instead we are able to test it. In addi- 
tion to this, Efstathiou assumes that cooling in the interior of bub- 
bles inflated by SNe is negligible and therefore SNe remnants can 
only contribute to the hot phase of the ISM. In our model we allow 
the interior of bubbles to cool down, which is a key process to fol- 
low, as in most of the cases cooling is efficient and bubbles enter a 
radiative phase rather quickly. 

We find that our predi c ted ou tflow rates are similar to those 
predicted by iHopkins et alj ( 1201 2h using simulations that resolve 
scales just below the size of GMCs and model SNe feedback by 
injecting thermal energy stochastically into neighbouring particles. 
However, their outflow rates correspond to the sum of several pro- 
cesses, such as photoevaporation and radiation pressure, and are 
not exclusively SNe driven outflows. They argue that in dense en- 
vironments, radiation pressure dominates the overall outflow rate. 
In those environments our scheme predicts a larger contribution to 
the outflow rate from SNe than that predicted by Hopkins et al. 
Nonetheless, note that we indirectly assume that photoionisation 
takes place due to our assumption of SNe driving bubbles which 
expand against the warm medium instead of the dense gas from 
which stars form. 



5 TOWARDS A NEW PARAMETRISATION OF THE 
OUTFLOW RATE 

One of the main aims of this paper is to establish if the results of 
our dynamical model of SNe feedback can be reproduced using a 
simple parametrisation cast in terms of global galaxy properties. In 
this section we use our dynamical model of SNe feedback embed- 
ded in GALFORM to assess parametrisations of the mass loading 
used in the literature (§ 15. Il l and search for an improved way of 
reproducing the mass loading factor (ij |5.2t . 



5.1 Dependence of the outflow rate on circular velocity 

As discussed in the Introduction, a widely used approach in galaxy 
formation models is to parametrise the mass loading of the outflow 
solely in terms of the circular velocity, Hcirc, which is considered as 
a proxy for the depth of the potential well of the galaxy. Scalings of 
P with circular velocity can be motivated by invoking momentum- 
conserving (/3 cx v~^l^) or energy-conserving (/3 oc u^^^,) winds, 
or the power-law index can be treated as a free parameter, as in 
GALFORM and most other semi-analytic models. Our model has the 
power to test such assumptions by directly comparing the j3 calcu- 
lated for a given timestep with the circular velocity of the galaxy. 

Parametrisations of SNe feedback that include a direct scaling 
with the circular velocity of the galaxy can be grouped into two: 
those assuming a single scaling relation for both the outflow rate 
from the galaxy and from the halo, and those which separate them 
into two different mass loading factors, /3ism for the mass loading 
of the galaxy and /3haio for that of the halo. G ALF ORM is an e xample 
of the first type (see also iLagos et aflboOSl a ndlCook et al 



In the second t ype, we find the mod e ls of e .g. lCroton et al.l 
[Monaco et al] OOOTh, iMaccio et alj j2010h and IGuo etal] 




Croton et al. ( 2006h assume that the outflow rate from 



For instance, 

the galaxy sc ales linearly with th e ins tantaneous SFR, and adopt 
/3isM = 3.5. lMacci6etal.l ( l2010h and lOuo et J] j201 ih modified 
the form of /3ism so that it makes a transition from a constant value 
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function of the circular velocity of the disk for quiescent galaxies with 
Mi, > 10* h~^MQ in the model with the standard choice of parameters 
(see Table[T). The relation is shown for different redshift ranges, as labelled. 
Solid lines and errorbars indicate the median and 10 and 90 percentile 
ranges of the relations. We also show the parametrizations used in a range of 
semi-analytic models, corresponding to (i) Baugh et al. (2005; dotted Hne), 
(ii) Dutton et al. (2010; dot-dashed line), (iii) Bower et al. (2006; solid line) 
and (iv) Guo et al. (2011; dashed line) (see text for details of the models). 
Bottom panel: The /3 — f circ relation in the model with the standard choice 
of ISM parameters for starburst galaxies with M* > 10* /i~^Mq at dif- 
ferent redshifts. In this case the circular velocity corresponds to that of the 
bulge. Lines and colours have the same meaning as in the top panel. 



dashed line in Fig. 1131). In the Dutton et al. model, the normalisation 
velocity is calculated from the momentum injected by a single SN 
that ends up in the outflow, which is 3.2 x 10* Mq kms^^ for a 
Kennicutt IMF. 

(iii) P = (t;circ/485kms~^)~^-^ from lSower et al.l ( |2006|) (solid 
line in Fig. 1131). 

(iv) /3 = 6.5 [0.5 + (ucirc/70kms-^)-^-^] from lGuo et ij] ( l201ll) 
(dashed line in in Fig. I13t . which gives a SNe driven wind with a 
high mass loading even in galaxies with veiy high circular veloci- 
ties, e.g. corresponding to those at the centre of clusters. 

There are three key conclusions that can be drawn from 
Figim (i) 3 single power-law fit cannot describe the dependence 
of P on Hcirc, (ii) there are large variations in the normalisation, 
but also in the slope of the /3-iicirc relation with redshift, and (iii) 
starbursts and quiescent galaxies follow different relations. 

Regarding the shape of the P-Vdrc relation, the top panel 
of Fig[T3] shows that our dynamical calculations display a trend 
of P decreasing with increasing «circ for galaxies with Hdrc 
SOkms^^. Below «circ ~ 80kms~^, the predicted mass load- 
ing shows a flattening or even a turnover followed by a positive 
P-Vcirc relation. The parametrisations used in the literature for the 
relation between P and Ucirc, are a poor description of the relation 
obtained from our physical model, which does not display a simple 
pow er-law behaviour w hen plotted in this way. 

iFont et al] l l201lh discuss a phenomenological model with a 
saturation of the SNe feedback, which was invoked to reproduce 
the observed LF and metallicity of the Milky Way's satellites. Font 
et al. set a ceiling P = 620 for Ucirc < 65kms^^ to obtain a 
good match to the properties of the Milky Way's satellites. Our 
dynamical model of SNe feedback predicts a qualitatively similar 
behaviour to the saturated feedback scheme of Font et al. The peak 
value of /3 at z = is similar to the saturation value proposed 
by Font et al. However, we find that the peak value of the mass 
loading and the circular velocity at the peak occurs change with 
redshift. We also find that saturation velocity varies with the param- 
eters adopted to describe the ISM and molecular clouds, spanning 
the range Hdrc.sat ~ 70 — 100 km . In our model the saturation 
velocity has no direct connection to the ratio between SNe energy 
and halo potential. 

The redshift variation of the mass loading of the wind 
can be quantified by fitting a power law of the form P — 
(I'circ/I^iot)""'"' to quiescent galaxies at different redshifts 
(top panel Fig|13b. For circular velocities in the range Vdic > 
80km s~^, the dependence of cihot and Vhot on redshift is given 
by 



in high circular velocity galaxies to a form in which Pism increases 
as the circular velocity of the galaxies decreases, in order to better 
reproduce the number density of low-mass galaxies (see (iv) is the 
list below). In our model we calculate Pism and compare it with 
the parametrisation from 4 of the previous models. 

Fig[T3]shows the /3 predicted by the dynamical SNe feedback 
model after implementing it in the full galaxy formation simulation, 
plotted as a function of circular velocity for quiescent (top panel) 
and starburst galaxies (bottom panel). The model shown in Fig|13l 
corresponds to the standard choice of model parameters (see Ta- 
ble[T]). We overplot for comparison the following parametrisations 
for the mass loading from the literature; 

(i) P = («circ/300kms-^)-^ from iBaugh eTal] ( |2005|) (dotted 
line in Fig. 1 131). 

(ii) P = (ucirc/300kms"^)"^ from button et al.l ( 1201 Ol) (dot- 



Qhot = 2.7 + 21og(l + 2), (55) 
l^hot = 425kms"' (l + z)~°-^ (56) 

For galaxies with Ucirc/kms~^ < 80 and for starbursts, the de- 
pendence of Qfhot and Vhot on redshift is more complicated and 
cannot be described by simple power-law fits. This behaviour illus- 
trates that the mass loading of the outflow does not have a natural 
dependence on circular velocity. 

When focusing on starburst galaxies only, we find that the de- 
pendence of P on Ucirc changes dramatically (see bottom panel of 
Fig. 1131 ). This is due to the very different conditions in the ISM 
in starbursts compared to quiescent galaxies, with higher gas sur- 
face densities for a given Ucirc- The turnover obtained for quiescent 
galaxies at Vdrc ~ 80kms~^ is also present in starburst galax- 
ies at z < 2. We find that the differences between quiescent and 
starburst galaxies and the turnover at Vdic ~ 80kms~^ can be 
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explained in terms of the more fundamental relation between /3 
and the gas scaleheight, hg. For the latter case, both quiescent and 
starburst galaxies follow nearly the same relation (see top panel 
of Fig. [8}. This explains the nature of the P-Vdic relation: there 
is a correlation between t;circ and hg, for quiescent galaxies with 
'^circ > 80kms~^, but this is not present at lower Vdrc or in star- 
burst galaxies. 



5.2 A new parametrisation of the mass outflow rate 

We analyse the dependence of j3 on various properties of the disk in 
order to find the most natural combination of parameters to describe 
the mass loading. This new way of describing /3 can therefore be 
used in semi-analytic galaxy formation models and simulations. 

Fig. [14] shows the mass loading factor, /3, as a function of (i) 
Egas, (ii) Pgas, (iii) Egas + E* and (iv) ftg, for the standard set of 
parameters for GMCs and the diffuse medium (see Table [T). Note 
that the third of these quantities can be written in terms of the sur- 



Standard model 



face density of gas and the gas fraction Eg 



^gas 

//gas. All 



quantities above are evaluated at the half-mass radius of the disk or 
the bulge, rso (see Appendix iBl for the definition of the profiles), 
and the predictions are shown for all galaxies, quiescent and SB, in 
different redshift ranges. We decide to study the relation between 
P and these quantities due to the correlation we find between the 
mass of a single bubble at the point of break-out from the disk and 
the local properties Pgaa, Egas, Egas + E* and /ig (see Fig.|5]l. We 
also show the resulting relation between /? and the quantity plotted 
on the a;-axis if we use the old mass loading parametrisation (see 
point (iii) in list of ij IS.lb . 

We find that our results can be approximately described by the 
following fits 



/3 = 
/3 = 



Egas(r5o) 



1.6 X 103 Mq pc-2 

/ \ 1 -0.5 

Pgas (rso ) 



14 Mq pc- 



,(r5o) + E4r5o) 



2.6 X 103 Mq pC-2 

hgirsoy ^'^ 



spc 



(57) 
(58) 
(59) 
(60) 



We quantify how good the correlation is by using two statis- 
tics, the Pearson correlation coefficient, R, and an estimate of the 
dispersion around the median, (Tm. For each a;-axis bin we calculate 
a dispersion, (Jx, corresponding to the ratio between the sum of the 
square of the deviations around the median in the y-axis and the 
number of objects in the bin. We then calculate am, which corre- 
sponds to the square root of the median value of the distribution of 
iTx. We calculate am in the log-log plane, in units of dex. Note that 
R and am are independent statistics which can be used to assess 
how good the correlation is between two quantities. The values for 
both quantities for galaxies at z < 0.1 are written in each panel of 
Fig. [11] 

In terms of the Pearson correlation factor, R, and the dis- 
persion, am (shown in Fig. I14t . the properties that best describe 
/3 are Egas + ^* and /ig. Fig. [14] shows that the normalisation 
and power-law index of the above relations vary with redshift, 
with high-redshift galaxies following a steeper relation than low- 
redshift galaxies. This trend can be understood as being due to high- 
redshift galaxies having larger gas fractions compared to lower- 
redshift galaxies. Galaxies with a high gas fraction typically have 



"a 

0) 




-4 -2 2 4 

log(p fit) 

Figure 15. The predicted mass loading from the full model (y-axis) plotted 
against the fit given by Eg. 1621 expressed in ternis of the gas scaleheight and 
gas fraction, for galaxies at different redshifts, as labelled. 



a molecule-dominated ISM, and these are predicted to follow a 
steeper relation between /3 and hg than those with an atomic- 
dominated ISM, which are typically gas poor (see § 14.1.11 for an 
analytic derivation of such a trend). We find that the redshift trend 
can be removed by adding an extra dependence on the gas fraction 
to the expressions for /3, 



Egas(r5o) 



1600 Mq pc- 



■ -0.6 


/gas 







0.12 







1.1 


/gas 


0.4 


15 pc 




0.02 





(61) 
(62) 



which both have Pearson correlation factor of i? « 0.97 and a dis- 
persion am ~ 0.3 dex for galaxies at z < 0.1. This is shown in 
Fig. [15] where the fit of Eq.|62]is compared with the directly calcu- 
lated p. Most of the redshift evolution seen in Fig.[T4]is removed. 

Eqs. [61] and [62] are also useful to characterise the mass 
loading /3 obtained in the model when varying the parameters 
used in the ISM modelling (Table [TJ. This is shown in Fig. |16| 
in which the power-law indices and normalisations for the re- 
lations, defined as /3 = (Egas/Eo)''« (/gas//o,E)''f « and = 
(/ig/fto)'''' (,fgas/./o,h)'''''\ are shown for 3 different choices of 
ISM model parameters. The model using fi — 1.1 corresponds 
to the weakest feedback model and that with usf = 0.3 Gyr^^ 
to the strongest feedback model. The three choices of model pa- 
rameters produce very little variation in the power-law indices of 
the above relations (top panel of Fig.ll6t. Variations are observed 
in the normalisations of the relations and represent different feed- 
back strengths (bottom panel of Fig. I16t . This means that if we 
were to include the parametric form given by Eqs.[6T]and[62]in the 
semi-analytic model, we would need to vary the zero-point of these 
relations to reproduce the results for different parameters for the 
diffuse ISM and GMCs. Eqs.l61landl62ldescribe our results for the 
mass loading /3 in galaxies at any redshift, within the range tested 
(i.e 2 < 10 and + Mgas.iSM > 10® /i"^ Mq) with very little 
dependence on redshift or stellar mass. 

The old parametrisation (shown by the dashed lines in Fig.ll4t 
results in a trend of /? decreasing with the properties plotted on the 
X-axis, given the correlation already discussed between Vdic and 
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Figure 14. The global mass loading factor, f) = A/pjcct/^, as a function of gas surface density, Sgas (top left-hand panel), gas density, pgas (top right- 
hand panel), gas plus stellar surface density, Sgas + (bottom left-hand panel) and the gas scale height, hg ( bottom right-hand panel), for galaxies with 
> 10® h~^M(T,. All quantities plotted on the x-axis are calculated at the half-mass radius of the disk in the case of quiescent SF, or the bulge in the 
case of starbursts. The relations are shown for different redshift ranges, as labelled, and correspond to the predictions of the model with the standard choice of 
parameters (listed in Table[T). Solid lines and errorbars indicate the median and 10 and 90 percentile ranges of the relations. For reference, the values of the 
Pearson correlation coefficient, R, and the dispersion aro und the median, am I dcx, calculated for galaxies at z < 0.1 in the new model are written on each 
panel. We also show the results obtained when using the iBower et al.l (2006) choice for the outflow rate, = (tJcirc /485 km s ~ ^ ) ~ ' ^ , for galaxies at 
2 < 1 and 6 < 2 < 8 (dashed lines) in each panel. The horizontal shading represents the 10 and 90 percentile ranges of the relations using the Bower et al. 
parametrisation. 



these variables. However, /3oid differs from the mass loading /3 for 
galaxies with low surface densities of gas by up to a factor of ~ 5 
in either direction, and overestimates /3 at the high surface den- 
sity regime by up to a factor of ~ 100, depending on the redshift. 
In Fig. [14] /3oid varies with redshift much more strongly than the 
new parametrisations, and therefore overestimates the SNe feed- 
back in high-redshift galaxies. This reflects the importance of the 
analysis performed in this paper and the need for a revision of such 
parametrisations. The largest differences between the predicted /? 
and /3oid are obtained at high-redshifts. 

The difference between SBs and quiescent galaxies apparent 
in the /3— Vcirc plane in Fig.[T3]is greatly reduced in the /3— /ig plane 
(see the top panel of Fig.[8l(. This is because SB galaxies of a given 
iicirc have much higher densities in stars and gas than their quies- 
cent counterparts. Although the relation is noisier due to the lower 
numbers of SBs in the model output compared to quiescent galax- 
ies, the j3 — hg relation is very similar in slope and normalisation 
to that for quiescent galaxies. This suggests that the dependence 



of mass loading is fundamental and captures the relevant physics 
determining /3. 



6 THE IMPACT OF THE NEW OUTFLOW MASS 
LOADING ON GALAXY FORMATION 

In this section we consider the impact of our dynamical model of 
SNe feedback on galaxy properties and compare with the predic- 
tions of the model which uses the old parametrisation. We first es- 
timate the error associated with using the parametric form defined 
in Eq. [62] instead of performing the full calculation caiTied out in 
this paper. Second, we analyse the net effect of our dynamical mod- 
elling on galaxy properties by focusing on two statistical properties 
of galaxies: (i) the evolution of the LF in the K- and V-bands, and 
(ii) the evolution of the global SFR density. An analysis of a com- 
plete set of galaxy properties will be presented in a future paper 
(Lagos, Lacey & Baugh, in prep.). 
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Figure 16. Top panel: Power-law slope of the relations between /3, /ig, 
Sgas and /gas, quantified as /3 = [h^/ho)''^ (/gas//g,h)''*'''' and /3 = 
(Sgas/So)''s (/gas//g,E)'''''' ■ Lines are as labelled in the panel. The re- 
sults for the model with the standard set of parameters and those predicting 
the highest and lowest /3 are shown as symbols. The parameters of the fit 
correspond to fitting the relations above in subsamples of galaxies at z < 8 
with different gas fractions. Bottom panel: Normalizations of the relations 
above for the three models of the top panel. Lines are as labelled in the 
panel. The plot shows that the power-law slopes are not affected by changes 
in the parameters describing the ISM and SF but that only the normalisa- 
tions of the relations change. 



We ran the full dynamical model in which /? is calculated self- 
consistently, and compare with the model using the prescription 
from Eq. |62]to calculate /3, under the simplifying assumption of 

= /3. We compared the luminosity functions predicted by both 
procedures in the bands 900 — 1200A, bj, V, K and 8^m. At 
z = 0, the largest differences are obtained in the far-UV band, 
but are at most ~ 25%. The other bands show differences in the 
range 5 — 20%. However, at z = 6 these differences can be as 
large as 80%. The reason for the larger differences at high redshifts 
is that we currently do not allow for variations in the parametri- 
sation of fi^ with respect to /?, like those shown in Fig. [TT] Such 
variations have only a minor effect at z = 0, but they have an af- 
fect in z > 4 galaxies, where larger differences between /? and 

are predicted by the dynamical model. The main drivers of the 
differences seen in the luminosity functions are differences in the 
cold gas mass and mass in metals in the ISM. The stellar mass 
and hot gas mass functions are similar to within < 40% at red- 
shifts z = — 6. In the redshift range shown in Figs.ll7landll8l 
variations between the self-consistent calculation and the calcula- 
tion using the /3 parametrisation are not significant. We calculate 
the best parametrisations using the form of Eq. |62]for the different 
ISM parameter choices and present in Table |3]the results for four 
choices of parameters spanning the full range of feedback strength. 
We find that using the prescription for P given in Eq.|62]gives reli- 
able results that closely follow the behaviour of the full dynamical 
model at z < 4, but significantly speeds up the calculation. 

In order to analyse the effect of the new dynami cal model of 
SNe feedback on galaxy properties, we focus on the iLagos et al.l 



Table 3. Models shown in Figs. [Tt] [Hand [19] The first row gives the old 
parametrisation used to describe the outflow. The next four rows show alter- 
native models using the new /3 parametrisation of Eq. |62| Each parametri- 
sation represents different parameter choices for the full SNe feedback dy- 
namical model, which is indicated in the parenthesis. The parametrisation 
used for each model is shown in the second column. 



Model 



/3 parametrisation 



V485 kms-l J 



. 0, 



Lagosl2.01dBeta V 485kn.s-i 

Lagosl2.WeakSN(/, = 1.1) [23vcJ \ o.S 

Lagos 12.InterSNa (riif;, qmc = 0-03 Gyr) ( 
Lagos 12.InterSNb (Std.) 

Lagosl2.StrongSN(i.sP = 0.3 Gyr"!) (;^)''' (w) 



0.4 
0.4 



17 pc/ V 0.1 J 

15 pc/ V0.02/ 

1.1 / f \ 0A 



( I2OI2I) model and vary the SNe feedback prescription. We compare 
the four alternative models listed in Table [3] 

Fig. [17] shows the A'-band LF at various redshifts for the 5 
models listed in Table |3] We remind the reader we are not trying 
to fit observations here, but rather we are trying to see the effect 
the modelling of feedback has on galaxy properties starting from a 
model which uses a completely different way of calculating /?. The 
most interesting feature in Fig.[T7]is that all the models that use the 
new feedback model developed in this paper give a shallower faint- 
end slope at 2; < 2.5, regardless of the ISM model parameters, but 
produce a higher overall normalisation for the LF. The model with 
the strongest feedback (Lagosl2.StrongSN) shows a faint end that 
is similar to the original model. There is a trend of a shallower faint 
end with weaker SN feedback models, although this trend changes 
with band and redshift. It is also clear that the models predict very 
weak evolution of the slope of the faint end. The shallower faint 
end slope predicted by our new feedback scheme suggests that the 
problem of the predicted steep faint end of the LF and low-mass 
end of the stellar mass function could be largely overcome by using 
the new parametrisation of the mass loading (Eq.l62). The physical 
reason behind the shallower faint end slopes obtained by using the 
new /? parametrisation is that faint galaxies typically have large /ig 
and therefore can reach very large values of /3. These faint galaxies 
do not necessarily correspond to those with the smallest Ucirc, and 
therefore in these galaxies, the new parametrisation drives larger 
P than that obtained with the Vdic parametrisation. We show in 
Fig.[T4]that in the regime of large /ig, the parametrisation used by 
the original model (see Table [3} gives a lower /? compared to the 
new parametrisation. 

The bright-end of the if-band LF predicted by the models us- 
ing the new feedback prescription is higher in all the cases com- 
pared to the original model. This is due to the lower /? predicted 
by the dynamical SN feedback model compared to the parametri- 
sation adopted in the Lagosl2.01dBeta model. This, in addition to 
the unchanged gas reincorporation timescale, leads to more bright 
galaxies. In paper II we will model the expansion of bubbles in the 
halo to remove this process as a free parameter. We will analyse in 
more detail the effect of SN feedback on the bright end of the LF. 

Fig. [Tsjis equivalent to Fig.[T7]but shows the V^-band LF for 
z > 0.5. The behaviour of the models in this band is broadly the 
same as in the near-IR: the new feedback scheme, regardless of the 
strength of the SN feedback, predicts a shallower faint end of the LF 
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Figure 17. Rest-frame _ft'-band galaxy luminosity function for the Lagos 12 
model with the old and the new SNe prescriptions (see Table |3], at vari- 
ous redshifts, as labelled. Observational res ults from Pozzetti et al. 1 2003), 
iDrorv et al.|j2004 . ls^acco et al.N2006l) and lCaputi et al.. 12006.) are shown 
as grey symbols, identified by the key in the two top panels. Note that the 
models have not been retuned to fit the observed LF. 



up to 2 « 1.5. However, above that redshift, the strength of the SN 
feedback plays an important role in determining whether the faint 
end is shallower or steeper than predicted by the original model. 
The slope of the faint end in the V-band LF varies more strongly 
with redshift and in a complex way compared to the variations seen 
in the iC'-band LF. 

Interestingly, the different SNe feedback models of Table |3] 
converge to similar LFs in both the K and V bands at 2; > 3 but 
evolve differently towards z = 0. This is because these models pre- 
dict galaxies with different star formation histories. Fig.|19|shows 
the global SFR density evolution predicted by each of the models 
of Table[3l The models using the new SN feedback scheme predict 
that the global SFR peaks at slightly lower redshifts compared to 
the original model, with weaker SN feedback producing a lower 
redshift for the peak. Note that even the model with the strongest 
SNe feedback produces larger SFR densities at z « 2 — 4 compared 
to the model using the old jS parametrisation. Compared to obser- 
vations, the model with the strongest SN feedback predicts SFR 
densities that are too low, while the weakest SN feedback give SFR 
densities that are too high. It is interesting to note that the model 
with the strongest SN feedback results in the largest decline in the 
global SFR per unit volume, dropping by a factor of ~ 30 from 
the peak to the present day. A key physical process to analyse be- 
fore ruling out any of these models is the reincorporation timescale 



V-band 
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Figure 18. Rest-frame V-band galaxy luminosity function for the Lagos 1 2 
model with the old and the new SNe prescriptions (see Table[3), at various 
redshifts, as labelled. Observational results from iMarchesini et alj j2012l) 
are shown as grey symbols. Note that the models have not been retuned to 
fit the observed LF. 



of the gas after outflowing from the ISM into the hot gas reservoir 
of the halo. Also, other galaxy formation parameters may have to 
be reset, since these were based on the old outflow model. In any 
case, the fact that the use of the new /? parametrisation predicts a 
shallower LF of galaxies points to the need to revise the physics 
included in galaxy formation models and simulations. 



7 DISCUSSION AND CONCLUSIONS 

We have presented a dynamical model of SNe feedback which 
tracks the evolution of bubbles inflated by SNe into the ISM of 
galaxies. Our model includes a range of processes which can affect 
the expansion of bubbles: gravity, radiative energy losses, exter- 
nal pressure from the diffuse medium and temporal changes in the 
ambient gas. Bubbles inflated by SNe are evolved from the adia- 
batic to the radiative phases until the point of break-out from the 
galaxy disk or bulge, or confinement in a multi-phase ISM. The 
multi-phase model of the ISM includes a diffuse, atomic phase, a 
dense, molecular phase and a hot, low density phase. The latter 
corresponds to the interior of bubbles. The metal enrichment of the 
ISM and halo due to SNe takes place through bubbles. The loca- 
tion of star-forming regions, or GMCs, which give rise to bubbles 
is connected to the radial distribution of molecular gas, which al- 
lows us to study both the global outflow rate and the radial profile 
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Figure 19. The evolution of the cosmic star formation rate per unit volume 
for the Lagos 12 model with the old and the new SNe prescriptions which 
give rise to different strengths of SN feedback (see Table [5], as labelled. 
The observational estimates of Karim et al. (2011; asterisks) an d the data 
compi lation of Hopkins et al. (2004; diamonds) are al so shown. iHopkinj 
<2004 assumes a Salpeter IMF and lKarimerS] iioi lh a Chabrier IMF. 
Therefore, SFRs have been scaled to a Kennicutt IMF (scaled down by a 
factor of 2 in the Salpeter case and down by a factor 1.12 in the Chabrier 
case). 



of galactic outflows. The aims of this work are (i) to test the impor- 
tance of each of the physical processes included in the expansion 
of bubbles and to explore the parameter space of the modelling of 
GMCs and the ISM, (ii) to determine which combinations of galaxy 
properties the outflow rate best correlates with and (iii) to improve 
upon widely used parametric forms for the outflow rate used in the 
literature. 

To help us assess these points, we embed our calculations in 
the GALFORM semi-analytic model, which follows the formation 
and evolution of galaxies in the framework of hierarchical structure 
formation. We take advantage of the two-phase me diu m description 
introduced into GALFORM bv lLagos et al J ( |20 1 id) and lLagos et al.l 
( 12011 j) . to trace star formation and star forming regions using 
the cold molecular component of the ISM, while allowing bub- 
bles to sweep up gas only from the diffuse neutral atomic com- 
ponent. In the Lagos et al. model, the molecular-to-atomic mass 
ratio is calculated from the radial profile of the hydrostatic pres- 
sure, and the SFR is calculated from the molecular g as radial pro- 
file (e.g. iBlitz & Roso lo wskv ,20o3 ; iLerov et al.ll2008.) . The semi- 
analytic model provides the initial conditions needed by the dy- 
namical model of SN feedback: the stellar and dark matter contents, 
the surface density of atomic and molecular gas, the gas metallic- 
ity and the scalelength of each mass component. This modelling 
allows us to study the relation between the rate at which mass es- 
capes from the galaxy disk or bulge (outflow rate) and the proper- 
ties of the disk, bulge and halo, over a wide dynamic range. Pre- 
vious work has focused on hydrodynamical simulations covering 
a narro w dynamic ranee, which has been chosen somewhat arbi- 
trarily dHopkins et al.l 20121 ; ICreasev et al.ll2013h , or which have 
adopted Sedov analytic solutions for the evolution of bubbles (e.g. 
lEfstathiou"2000l; lMonacoll2004 p). One of our goals is to comple- 
ment and extend this work by using a more general SNe feedback 
model and the galaxy population and star formation histories pro- 
duced by the semi-analytic model. 



We summarise our main conclusions below: 

(i) We find that the mass loading of the outflow, decreases 
with increasing gas surface density and increases with increasing 
gas scaleheight. On the other hand, the outflow velocity increases 
with increasing gas surface density and decreases with increasing 
gas scaleheight. These trends are seen in both the global and local 
mass loading and velocity of the wind. 

(ii) We find that the multi-phase ISM treatment included in 
our model is essential for reproducing the observed outflow rates 
of galaxies. When fixing the diffuse-to-cloud mass ratio instead of 
calculating it from the hydrostatic pressure, we find variations in 
the predicted mass loading /3 of up to 2 orders of magnitude in the 
highest gas density regimes. This emphasizes the importance of the 
multi-phase ISM included in our modelling. By adopting different, 
but still plausible parameters in the modelling of GMCs and the 
diffuse medium, we find variations in /3 of a factor up to ~ 3 and 
in ^outflow of a factor up to ~ 1.7 in either direction. We also find 
that by the time bubbles escape from the ISM, they are radiative in 
the majority of the cases. 

(iii) When comparing our predicted outflow rates and ve- 
locities with those inferred from observations (e.g. lMartiii|[l999l ; 
iBouche et alj20l3) . we find good agreement. We also find that our 
predictions are similar to those from the non-co smo logical hydro- 
dynam ical simulations of Hop kins et al.l j2012h and ICreasev et al.l 
1 I2OI3I) , in the regimes they were able to probe. Our work therefore 
confirms the finding that the surface density of gas is an important 
quantity in determining the mass loading of the outflow. 

(iv) The widely used parametric forms describing SNe feed- 
back and relating the mass loading /3 to only the circular velocity 
of the galaxy do not capture the physics setting the outflow rates 
from galaxies. For instance, we find that the trend of /3 decreasing 
with Ucirc is only valid for galaxies with Udrc ^ SOkms^^. Be- 
low this threshold, j3 flattens or decreases with decreasing fcirc. We 
also find that the relation between jS and Vdrc changes substantially 
with redshift. We find that tighter relations are those between j3 and 
the gas scaleheight and gas fraction, /3 oc [/ig(r5o)]^'^ [/gas]*^'^, 
and between /? and the surface density of gas and the gas frac- 
tion, P oc [Egas(r5o)]~°'®[/gas]° **. Changing the parameters in 
the model of GMCs and the diffuse medium can change the nor- 
malisation of these relations, but does not alter the power-law in- 
dex. We find that starburst and quiescent galaxies follow similar 
relations, with starbursts being slightly offset to lower (3 compared 
to quiescent galaxies. The outflow velocities can also vary between 
starbursts and quiescent galaxies depending on the adopted star for- 
mation law. A more rapid conversion from gas to stars drives larger 
velocities due to the higher energy and momentum injection rate 
from SNe. 

(v) We study the effect of the dynamical model of SN feedback 
developed here on galaxy properties and test the inclusion of the 
new parametrisation of /? (see (iv) above). We find that the faint end 
of the near-infrared LP becomes shallower in the model using the 
new feedback scheme compared to the old model. We find that this 
shallowing of the faint end takes place regardless of the parameters 
assumed to describe the diffuse ISM and GMCs, with a trend of 
weaker SN feedback predicting a shallower faint end of the LP. 

Our model is subject to simplifications required to model the 
evolution of bubbles in the ISM of galaxies. A critical simplifica- 
tion we make is to fix the GMC mass. A more sophisticated ap- 
proach would be to include a distribution of GMC masses and their 
spatial distribution following a theoreti cal estimate of the spatial 
clustering of GMCs of different masses ( lHopkinsll2012l) . However, 
such a description also requires more detailed information about 
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the ISM. Instead, we test our predictions by varying tlie adopted 
GMC mass in the range allowed by observations (see Table[T), and 
find variations in the normalisation of the mass loading described 
in (iv), but with little impact on the power-law indices. 

The agreement we find b etween ou r model and detailed hydro- 
dynamical simulations ( Hopkins et al.|[201 2: Creasev et al. 20131) 
suggests that we capture the relevant physics determining the rate 
at which mass escapes from the ISM of galaxies, despite the simpli- 
fications made in our modelling. The advantage of our calculations 
is that a much wider range of ISM conditions can be explored than 
is feasible in the more expensive hydrodynamical simulations. We 
have given predictions for the outflow rate for a very wide range 
in galaxy properties and cosmic epochs. The method developed in 
this paper also allows radial profiles of the outflow rate to be ob- 
tained. The new generation of integral field spectr oscopy instru- 
ment s, such as KMOS in the Very Large Telescope dSharples et al.l 
|2004) and the Sydney-AAO Multi-object Int egral field spectro- 
graph jCroom et al.ll2012l : iFogartv et al.ll2012l) will make the ob- 
servations of outflows routine in local and high-redshift galaxies, 
and will allow us to constrain our model observationally. 
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APPENDIX A: THE RECYCLE FRACTION AND YIELD 
OF DIFFERENT STELLAR POPULATIONS 

The number of SNe per solar mass of stars formed, t^sn, is calcu- 
lated from the IMF, (j){m) oc dA''(m)/dm, as, 

r?sN = / (j){m)dm, (Al) 

JmsN 

where msN = 8 Mq and mmax ~ 120 M©. Ror the iKennicutj 
( Il983l) IMR adopted here, ??sn = 9.4 x 10"^Mq ^ (in the case of 
a Salpeter IMR, t/sn = 7.3 x 10"^Mq ^). In §[2TT] we define the 
mass injection rate from SNe depending on the recycled fraction of 
massive stars, iisN. This recycled fraction also depends on the IMR 
as, 

RsN = / (m - mrcmii)<^(m) dm, (A2) 

where mrcm is the remnant mass. Similarly, we define the yield 
from SNe as 

PSN ~ / mi(m)0(m)dm, (A3) 

J msN 

where mi(m) is the mass of metals produced by stars of initial 
mass m. We use the stellar evolution models of Marigo (2001) and 
IPortinari et al.l ( fl998b to calculate the ejected mass from intermedi- 
ate and massive stars, respectively. Ror a Kennicutt IMR, we obtain 
RsN = 0.14 and PSN = 0.018. 



APPENDIX B: RADIAL PROFILES OF THE STELLAR 
AND DARK MATTER COMPONENTS AND THE 
MIDPLANE PRESSURE 

An important driver in the evolution of bubbles is the gravitational 
attraction exerted by the stellar and dark matter components. We 
describe here how we calculate the mass enclosed by a sphere of 
radius R located at a distance d from the centre of the galaxy. We 
perform our calculations of bubble evolution in shells in the disk, 
which defines d (see i]2.2.1t . 

The total stellar plus dark matter mass within a sphere of ra- 
dius R displaced by d from the centre of the galaxy coiTesponds 
to 

Mt{R, d) = AU{R, d) + Mom{R, d), (Bl) 
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where Mt{R,d) = M*_disk(-R, d) + M*,buigc(-R, rf) is the total 
stellar mass, M*,disk(^Z, d) and M«,buige(^?, d) represent the mass 
in the disk and the bulge, respectively, and A/dm(-R, d) the mass 
in DM, in all cases enclosed in R. We describe below how we 
calculate the variables of Eq. lBll 

Disk radial profile. We assume disks are well described by a radial 
exponential profile with a scale radius r^, which is related to th e 
half-mass radius as rscdisk = 1-67 rg teinnev & Tremain3l2008l) . 
We define the stellar surface density of the disk at a distance d from 
the centre as 



S*,disk(rf) = -T. 5-e 

III ri 



(B2) 



Here, ' disk is the total stellar mass in the disk. If the relevant 
sphere of radius _R is at a distance d from the centre, then the stellar 
mass in the midplane of the disk exerting the gravitational attraction 
on the bubble is approximately 



47rfl3 E^,disk(d) 
3 2/i* ' 



(B3) 



Here, /i* is the scale height of the stars, which we estimate from 
the scale radius o f the disk following the empirical results of 
lKregeletalJl2002l) , r^/h^, = 7.3. 

Bulge radial profile. The potential well of a galactic bulge, <&(r), 
can be well described by a Dehnen profile iDehneril 1993 ) with 
7b = 3/2 which closely resembles a ide VaucouleursI ( 1953 ) r^^* 
profile. 



$(r) 



GMl 



bulge 



2 -7b 



r + ro 



(B4) 



where ro is the scale radius and A/,[: ^uige is the total stellar mass 
of the bulge. The scale radius relates to the half-mass radius of the 
bulge, rso.b as 



r50,b = ro (2 



l/(3-7b) 



(B5) 



In this definition of potential well, the volume density profile 
of stars is. 



P*, bulge (r) = 



(3 - 7b) M^^bnlge 



(B6) 



47r rTb(r + ro)4-Tb ' 

Although the stars in the bulge follow a De Vacouleurs profile, the 
gas is assumed to be better characterised by an ex ponential pro- 
file, a s has been observed in early-type gala xies (e.g. lCrocker et al.) 
I2OIII ; iDavis etal.ll201ll : ISerraet al.ll2012h . This means that the 
same geometry adopted for the case of disks applies here: bubbles 
expand in a coordinate system displaced by d in the x-axis. How- 
ever, the difference with the case of the disk is that here the stellar 
profile has spherical symmetry. With this in mind, we approximate 
the stellar mass enclosed by a bubble of radius R displaced by d 
from the centre as. 



M,^, bulge (i?,d) 



P*,bulgc(rf), 



(B7) 



We use the equations above to calculate the Af* {R, d) that 
goes into Eqs. [Tl3][T8l20l and l24l26l 

Dark matter radial profile. He re we assume that D M halos are well 
described by a NFW profile dNavarro et al.lll997l) . We follow the 
description of ICole et al.l ( I2OO0I) . where halos contract in response 



to the presence of baryons. The galaxy disk, bulge and DM halo 
adjust to each other adiabatically. 

The volume mass density of DM is described in a NFW profile 

as 



pDM(r) = 



5c Pc 



(r/rs)(l+r/rs)2' 



(B8) 



where r^ is the DM scale radius, is the characteristic (dimen- 
sionless) density and pc is the critical density of the universe. As 
before, the mass enclosed within a sphere of radius R displaced by 
d from the centre of the potential well. 



AfDM(i?,d) 



4nR-^ 



PDM(d), 



(B9) 



assuming pomid) is approximately constant within the bubble. 

Note that Eqs. IB7I and IB9I are accurate in the regime where 
d/R ^ 1. In this paper we neglect the effect of tidal forces on 
bubbles, which arise from the asymmetric gravitational field, which 
distort their shape. This would affect the size of bubbles perpendic- 
ular to the gaseous disk and therefore the break-out of bubbles. 



Bl The midplane hydrostatic pressure of disk galaxies 

Under the assumptions of local isothermal stellar and gas layers, 
and (J* > CTgas, the midplane hydrosta tic pressure in di sks, Pext, 
can be approximated to within 10% by ( lElmegreen|[l993l) 



Pext(r) 



— GEgas 



(r) 



Sgas(r) + 



cr*(r) 



(BIO) 



where Egas and 



E,t are the surface densities of gas and stars at 
r, respectively, and erg and cr* give the vertical velocity disper- 
sion of the gas and stars. We assume a constant gas velocity dis- 
persion, CTd ~ lOkms"^ jLerov et al. 2008). By assuming that 
E* ^ Egas, o"*(r) = ^7rGh,tE*(r), where ft* is the stellar scale 
height. This approximation could break down for very high redshift 
galaxies, whose disks are gas dominated. In such cases, we assume 
a floor of (J* > CTg. 



APPENDIX C: CALCULATION OF SWEPT-UP, 
CONFINEMENT AND BREAK-OUT MASS RATES 

The contribution from bubbles to the rate of change of the mass 
and metallicity in the ISM and hot halo gas, depends on their 
evolution. In this appendix, we briefly describe how we calculate 
the overall contribution from bubbles in different evolutionary 
stages included in the set of Eqs. U6|45l 

The swept-up mass. Each galaxy has generations of bubbles 
whose evolution depends on the time they started their expansion 
and their spatial distribution in the galaxy. Each galaxy has its star 
formation history (SFH) sampled in a fine grid in time that goes 
down to the current time, tc. Each time interval, dt', in the SFH 
of a galaxy has associated a new generation of A^GMC,i,t' set of 
bubbles in the annulus i of the galaxy disk. Each of these bubbles 
have swept-up a mass msw(ri,t') from the diffuse medium and 
have a total mass mb(ri,t') at t' . The number of annuli used to 
solve the equations of bubbles expansion (§4.1) is N-^. The overall 
rate of swept-up mass is 
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Msw,ISM(tc) 




'cMCi.t' msw(n,t') (1 - Hr,h) 



(1 - At'. 



(CI) 



Here, Hr,h and Hv.a are step functions defined in terms of the 
radius of bubbles, R^, the gas scaleheight, /ig, the expansion speed 
of bubbles, Vs, and the velocity dispersion of the warm gas phase 
of the ISM, (Td, as, Hr,h = H[frhg{ri,t') - Rb{n,t')] and 
Hy^a ~ H[ad — Vs{ri,t')]. The quantities /ig, _Rb and depend 
on time and annulus. Eq. ICll implies that all bubbles contribute 
to the swept-up mass rate unless they have been confined or 
broken-out from the ISM in previous times. Bubbles at different 
evolutionary stages can coexist in an annulus. 

Confined bubbles. Confined bubbles contribute positively to 
A^g.iSM. The confinement of bubbles depends on whether the ex- 
pansion velocity of bubbles reach or exceed the velocity dispersion 
of the warm phase in the ISM, . The rate of mass transferred to 
the ISM by confinement is 



Break-out of bubble.^. The break-out of bubbles from the ISM 
contributes positively to the ISM gas due to the fraction of gas mass 
in the bubbles that stays in the ISM, (1 — /bo). The condition for 
break-out is that the radius of the bubbles reaches a factor /r of the 
gas scaleheight, i?b ^ frhg. The rate of break-out gas mass in the 
ISM is 



Afconf,ISM(ic) 




''GMC,i,t' mb(ri, t') Hy^a &t'{C2) 



A^bo,ISM(ic) 




GMC,i,t' mb(ri, t') Hr;h At' . (C3) 



